diff --git a/book/_config.yml b/book/_config.yml index f19fa9d..251f2d4 100644 --- a/book/_config.yml +++ b/book/_config.yml @@ -39,6 +39,7 @@ execute: - "**/geospatial-advanced.ipynb" - "cloud-computing/04-cloud-optimized-icesat2.ipynb" - "cloud-computing/atl08_parquet_files/atl08_parquet.ipynb" + - "machine-learning/photon_classifier.ipynb" allow_errors: false # Per-cell notebook execution limit (seconds) timeout: 300 diff --git a/book/_toc.yml b/book/_toc.yml index 8636460..9ef05d9 100644 --- a/book/_toc.yml +++ b/book/_toc.yml @@ -38,7 +38,7 @@ parts: - file: tutorials/cloud-computing/atl08_parquet_files/atl08_parquet options: - titlesonly: true - - file: tutorials/photon_classifier + - file: tutorials/machine-learning/photon_classifier.ipynb - caption: Projects chapters: - file: projects/index @@ -48,4 +48,3 @@ parts: - file: reference/bibliography - file: reference/IS2-resources - file: reference/questions - diff --git a/book/tutorials/index.md b/book/tutorials/index.md index 8f6d7a2..cfb40b4 100644 --- a/book/tutorials/index.md +++ b/book/tutorials/index.md @@ -10,4 +10,4 @@ Below you'll find a table keeping track of all tutorials presented at this event | [ICESat-2 Mission](./mission-overview/icesat-2-mission-overview.ipynb) | ICESat-2 Mission and Products | n/a | Not recorded | | [Cloud Computing](./cloud-computing/00-goals-and-outline.ipynb) | Cloud Computing Tutorial | n/a | Not recorded | | [Notebooks to Packages](./nb-to-package/index.md) | All about Python classes to packages | n/a | Not recorded | -| [ICESat-2 photon classification](./photon_classifier) | Machine Learning, PyTorch | ATL07 | Not recorded | +| [ICESat-2 photon classification](./machine-learning/photon_classifier.ipynb) | Machine Learning, PyTorch | ATL07 | Not recorded | diff --git a/book/tutorials/machine-learning/photon_classifier.ipynb b/book/tutorials/machine-learning/photon_classifier.ipynb new file mode 100644 index 0000000..849ba1f --- /dev/null +++ b/book/tutorials/machine-learning/photon_classifier.ipynb @@ -0,0 +1,14286 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "4ce203f6", + "metadata": {}, + "source": [ + "# Machine Learning with ICESat-2 data\n", + "\n", + "A machine learning pipeline from point clouds to photon classifications.\n", + "\n", + "Reimplementation of [Koo et al., 2023](https://doi.org/10.1016/j.rse.2023.113726),\n", + "based on code available at https://github.com/YoungHyunKoo/IS2_ML." + ] + }, + { + "cell_type": "markdown", + "id": "554bbfcf", + "metadata": { + "lines_to_next_cell": 2 + }, + "source": [ + "```{admonition} Learning Objectives\n", + "By the end of this tutorial, you should be able to:\n", + "- Convert ICESat-2 point cloud data into an analysis/ML-ready format\n", + "- Recognize the different levels of complexity of ML approaches and the\n", + " benefits/challenges of each\n", + "- Learn the potential of using ML for ICESat-2 photon classification\n", + "```\n", + "\n", + "![ICESat-2 ATL07 sea ice photon classification ML pipeline](https://github.com/user-attachments/assets/509dab2d-d25d-417f-97ff-fc966f656ddf)" + ] + }, + { + "cell_type": "markdown", + "id": "0e858676", + "metadata": {}, + "source": [ + "## Part 0: Setup\n", + "\n", + "Let's start by importing all the libraries needed for data access, processing and\n", + "visualization later. If you're running this on CryoCloud's default image without\n", + "Pytorch installed, uncomment and run the first cell before continuing." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "bf4f93e6", + "metadata": {}, + "outputs": [], + "source": [ + "# !mamba install -y pytorch" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "c475fb01", + "metadata": {}, + "outputs": [], + "source": [ + "import earthaccess\n", + "import geopandas as gpd\n", + "import h5py\n", + "import numpy as np\n", + "import pandas as pd\n", + "import pygmt\n", + "import pystac_client\n", + "import rioxarray # noqa: F401\n", + "import shapely.geometry\n", + "import stackstac\n", + "import torch\n", + "import tqdm" + ] + }, + { + "cell_type": "markdown", + "id": "6bc5754c", + "metadata": {}, + "source": [ + "## Part 1: Convert ICESat-2 data into ML-ready format\n", + "\n", + "Steps:\n", + "- Get ATL07 data using [earthaccess](https://earthaccess.readthedocs.io)\n", + "- Find coincident Sentinel-2 imagery by searching over a\n", + " [STAC API](https://pystac-client.readthedocs.io/en/v0.8.3/usage.html#itemsearch)\n", + "- Filter to only strong beams, and 6 key data variables + ancillary variables" + ] + }, + { + "cell_type": "markdown", + "id": "f6cdf9a9", + "metadata": {}, + "source": [ + "### Search for ATL07 data over study area\n", + "\n", + "In this sub-section, we'll set up a spatiotemporal query to look for ICESat-2 ATL07\n", + "sea ice data over the Ross Sea region around late October 2019.\n", + "\n", + "Ref: https://earthaccess.readthedocs.io/en/latest/quick-start/#get-data-in-3-steps" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "eddd2db0", + "metadata": {}, + "outputs": [ + { + "name": "stdin", + "output_type": "stream", + "text": [ + "Enter your Earthdata Login username: weiji14\n", + "Enter your Earthdata password: ········\n" + ] + } + ], + "source": [ + "# Authenticate using NASA EarthData login\n", + "auth = earthaccess.login()\n", + "# s3 = earthaccess.get_s3fs_session(daac=\"NSIDC\") # Start an AWS S3 session" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "39ef91fa", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "
\n", + " \n", + "
\n", + "
\n", + "
\n", + "
\n", + "
\n", + "

Data: ATL07-02_20191101183543_05480501_006_01.h5

\n", + "

Size: 42.6 MB

\n", + "

Cloud Hosted: True

\n", + "
\n", + "
\n", + " \"Data\"Data\n", + "
\n", + "
\n", + "
\n", + "
\n", + " " + ], + "text/plain": [ + "Collection: {'EntryTitle': 'ATLAS/ICESat-2 L3A Sea Ice Height V006'}\n", + "Spatial coverage: {'HorizontalSpatialDomain': {'Orbit': {'AscendingCrossing': 40.10013383465392, 'StartLatitude': -27.0, 'StartDirection': 'D', 'EndLatitude': -27.0, 'EndDirection': 'A'}}}\n", + "Temporal coverage: {'RangeDateTime': {'BeginningDateTime': '2019-11-01T19:39:35.790Z', 'EndingDateTime': '2019-11-01T19:54:06.084Z'}}\n", + "Size(MB): 42.600958824157715\n", + "Data: ['https://data.nsidc.earthdatacloud.nasa.gov/nsidc-cumulus-prod-protected/ATLAS/ATL07/006/2019/11/01/ATL07-02_20191101183543_05480501_006_01.h5']" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Set up spatiotemporal query for ATL07 sea ice product\n", + "granules = earthaccess.search_data(\n", + " short_name=\"ATL07\",\n", + " cloud_hosted=True,\n", + " bounding_box=(-180, -78, -140, -70), # xmin, ymin, xmax, ymax\n", + " temporal=(\"2019-10-31\", \"2019-11-01\"),\n", + " version=\"006\",\n", + ")\n", + "granules[-1] # visualize last data granule" + ] + }, + { + "cell_type": "markdown", + "id": "0108e5d4", + "metadata": {}, + "source": [ + "#### Find coincident Sentinel-2 imagery\n", + "\n", + "Let's also find some optical satellite images that were captured at around the same\n", + "time and location as the ICESat-2 ATL07 tracks. We will be using\n", + "[`pystac_client.Client.search`](https://pystac-client.readthedocs.io/en/v0.8.3/api.html#pystac_client.Client.search)\n", + "and doing the search with two steps:\n", + "\n", + "1. (Fast) search using date to find potential Sentinel-2/ICESat-2 pairs\n", + "2. (Slow) search using spatial intersection to ensure Sentinel-2 image overlaps with\n", + " ICESat-2 track." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "cc035a16", + "metadata": {}, + "outputs": [], + "source": [ + "# Connect to STAC API that hosts Sentinel-2 imagery on AWS us-west-2\n", + "catalog = pystac_client.Client.open(url=\"https://earth-search.aws.element84.com/v1\")" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "004d8e79", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + " 25%|██▌ | 2/8 [00:00<00:00, 12.12it/s]" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "b4fd119552fd40c089afb5b01a14bd9f", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "QUEUEING TASKS | : 0%| | 0/1 [00:00= 1:\n", + " # print(f\"Potential: {_item_len} Sentinel-2 x ATL07 matches!\")\n", + "\n", + " # 2nd check (spatial match) using centre-line track intersection\n", + " file_obj = earthaccess.open(granules=[granule])[0]\n", + " atl_file = h5py.File(name=file_obj, mode=\"r\")\n", + " linetrack = shapely.geometry.LineString(\n", + " coordinates=zip(\n", + " atl_file[\"gt2r/sea_ice_segments/longitude\"][:10000],\n", + " atl_file[\"gt2r/sea_ice_segments/latitude\"][:10000],\n", + " )\n", + " ).simplify(tolerance=10)\n", + " search2 = catalog.search(\n", + " collections=\"sentinel-2-l2a\",\n", + " intersects=linetrack,\n", + " datetime=f\"{start_time}/{end_time}\",\n", + " max_items=10,\n", + " )\n", + " item_collection = search2.item_collection()\n", + " if (item_len := len(item_collection)) >= 1:\n", + " print(\n", + " f\"Found: {item_len} Sentinel-2 items coincident with granule:\\n{granule}\"\n", + " )\n", + " break # uncomment this line if you want to find more matches" + ] + }, + { + "cell_type": "markdown", + "id": "7227bc62", + "metadata": {}, + "source": [ + "We should have found a match! In case you missed it, these are the two variables\n", + "pointing to the data we'll use later:\n", + "\n", + "- `granule` - ICESat-2 ATL07 sea ice point cloud data\n", + "- `item_collection` - Sentinel-2 optical satellite images" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8b6ffe42", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "369c2a5c", + "metadata": {}, + "source": [ + "### Filter to strong beams and required data variables\n", + "\n", + "Here, we'll open one ATL07 sea ice data file, and do some pre-processing." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "f06ac1ab", + "metadata": {}, + "outputs": [ + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "87e86f81abe346d3ba000a717472d911", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "QUEUEING TASKS | : 0%| | 0/1 [00:00" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "%%time\n", + "file_obj = earthaccess.open(granules=[granule])[0]\n", + "atl_file = h5py.File(name=file_obj, mode=\"r\")\n", + "atl_file.keys()" + ] + }, + { + "cell_type": "markdown", + "id": "8556669c", + "metadata": {}, + "source": [ + "Strong beams can be chosen based on the `sc_orient` variable.\n", + "\n", + "Ref: https://github.com/ICESAT-2HackWeek/strong-beams" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "ad7b0a60", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "['gt3r', 'gt2r', 'gt1r']" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# orientation - 0: backward, 1: forward, 2: transition\n", + "orient = atl_file[\"orbit_info\"][\"sc_orient\"][:]\n", + "if orient == 0:\n", + " strong_beams = [\"gt1l\", \"gt2l\", \"gt3l\"]\n", + "elif orient == 1:\n", + " strong_beams = [\"gt3r\", \"gt2r\", \"gt1r\"]\n", + "strong_beams" + ] + }, + { + "cell_type": "markdown", + "id": "ed79c372", + "metadata": {}, + "source": [ + "To keep things simple, we'll only read one beam today, but feel free to get all three\n", + "using a for-loop in your own project." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "479aaa97", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "gt3r\n", + "gt2r\n", + "gt1r\n" + ] + } + ], + "source": [ + "for beam in strong_beams:\n", + " print(beam)" + ] + }, + { + "cell_type": "markdown", + "id": "60abb52e", + "metadata": {}, + "source": [ + "Key data variables to use (for model training):\n", + " 1. `photon_rate`: photon rate\n", + " 2. `hist_w`: width of the photon height distribution\n", + " 3. `background_r_norm`: background photon rate\n", + " 4. `height_segment_height`: relative surface height\n", + " 5. `height_segment_n_pulse_seg`: number of laser pulses\n", + " 6. `hist_mean_median_h_diff` = `hist_mean_h` - `hist_median_h`: difference between\n", + " mean and median height\n", + "\n", + "Other data variables:\n", + "- `x_atc` - Along track distance from the equator\n", + "- `layer_flag` - Consolidated cloud flag { 0: 'likely_clear', 1: 'likely_cloudy' }\n", + "- `height_segment_ssh_flag` - Sea surface flag { 0: 'sea ice', 1: 'sea surface' }\n", + "\n", + "Data dictionary at:\n", + "https://nsidc.org/sites/default/files/documents/technical-reference/icesat2_atl07_data_dict_v006.pdf" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "ed923668", + "metadata": {}, + "outputs": [], + "source": [ + "gdf = gpd.GeoDataFrame(\n", + " data={\n", + " # Key data variables\n", + " \"photon_rate\": atl_file[f\"{beam}/sea_ice_segments/stats/photon_rate\"][:],\n", + " \"hist_w\": atl_file[f\"{beam}/sea_ice_segments/stats/hist_w\"][:],\n", + " \"background_r_norm\": atl_file[\n", + " f\"{beam}/sea_ice_segments/stats/background_r_norm\"\n", + " ][:],\n", + " \"height_segment_height\": atl_file[\n", + " f\"{beam}/sea_ice_segments/heights/height_segment_height\"\n", + " ][:],\n", + " \"height_segment_n_pulse_seg\": atl_file[\n", + " f\"{beam}/sea_ice_segments/heights/height_segment_n_pulse_seg\"\n", + " ][:],\n", + " \"hist_mean_h\": atl_file[f\"{beam}/sea_ice_segments/stats/hist_mean_h\"][:],\n", + " \"hist_median_h\": atl_file[f\"{beam}/sea_ice_segments/stats/hist_median_h\"][:],\n", + " # Other data variables\n", + " \"x_atc\": atl_file[f\"{beam}/sea_ice_segments/seg_dist_x\"][:],\n", + " \"layer_flag\": atl_file[f\"{beam}/sea_ice_segments/stats/layer_flag\"][:],\n", + " \"height_segment_ssh_flag\": atl_file[\n", + " f\"{beam}/sea_ice_segments/heights/height_segment_ssh_flag\"\n", + " ][:],\n", + " },\n", + " geometry=gpd.points_from_xy(\n", + " x=atl_file[f\"{beam}/sea_ice_segments/longitude\"][:],\n", + " y=atl_file[f\"{beam}/sea_ice_segments/latitude\"][:],\n", + " ),\n", + " crs=\"OGC:CRS84\",\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "818a77c4", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Total number of rows: 38246\n" + ] + } + ], + "source": [ + "# Pre-processing data\n", + "gdf = gdf[gdf.layer_flag == 0].reset_index(drop=True) # keep non-cloudy points only\n", + "gdf[\"hist_mean_median_h_diff\"] = gdf.hist_mean_h - gdf.hist_median_h\n", + "print(f\"Total number of rows: {len(gdf)}\")" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "9cf361ec", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
photon_ratehist_wbackground_r_normheight_segment_heightheight_segment_n_pulse_seghist_mean_hhist_median_hx_atclayer_flagheight_segment_ssh_flaggeometryhist_mean_median_h_diff
03.5128200.1289293470228.2500.08590338-53.941635-53.9435082.744329e+0700POINT (-166.17149 -66.1212)0.001873
14.1818180.1596233462216.0000.11759232-53.906929-53.9017372.744330e+0700POINT (-166.17152 -66.12129)-0.005192
24.1562500.1893043460210.2500.17018931-53.872688-53.8532182.744331e+0700POINT (-166.17154 -66.12139)-0.019470
34.5517240.2722343457500.0000.23654628-53.877270-53.8315662.744332e+0700POINT (-166.17156 -66.12148)-0.045704
44.6206900.2245593411361.2500.25243428-53.828976-53.7951892.744333e+0700POINT (-166.17159 -66.12156)-0.033787
.......................................
382417.7222220.1465301436920.2500.0754131718.94292818.9555403.303276e+0700POINT (18.94094 -63.78423)-0.012611
382427.9444450.1421201436920.2500.0482561718.92179518.9347383.303276e+0700POINT (18.94093 -63.78418)-0.012943
382437.5789480.1405041436920.2500.0337851818.89367518.9130253.303277e+0700POINT (18.94092 -63.78413)-0.019350
382446.7142860.1395101428857.8750.0269872018.86709018.9018973.303278e+0700POINT (18.9409 -63.78406)-0.034807
382456.8571430.1352621397495.8750.0093612018.85269518.8811913.303279e+0700POINT (18.94089 -63.784)-0.028496
\n", + "

38246 rows × 12 columns

\n", + "
" + ], + "text/plain": [ + " photon_rate hist_w background_r_norm height_segment_height \\\n", + "0 3.512820 0.128929 3470228.250 0.085903 \n", + "1 4.181818 0.159623 3462216.000 0.117592 \n", + "2 4.156250 0.189304 3460210.250 0.170189 \n", + "3 4.551724 0.272234 3457500.000 0.236546 \n", + "4 4.620690 0.224559 3411361.250 0.252434 \n", + "... ... ... ... ... \n", + "38241 7.722222 0.146530 1436920.250 0.075413 \n", + "38242 7.944445 0.142120 1436920.250 0.048256 \n", + "38243 7.578948 0.140504 1436920.250 0.033785 \n", + "38244 6.714286 0.139510 1428857.875 0.026987 \n", + "38245 6.857143 0.135262 1397495.875 0.009361 \n", + "\n", + " height_segment_n_pulse_seg hist_mean_h hist_median_h x_atc \\\n", + "0 38 -53.941635 -53.943508 2.744329e+07 \n", + "1 32 -53.906929 -53.901737 2.744330e+07 \n", + "2 31 -53.872688 -53.853218 2.744331e+07 \n", + "3 28 -53.877270 -53.831566 2.744332e+07 \n", + "4 28 -53.828976 -53.795189 2.744333e+07 \n", + "... ... ... ... ... \n", + "38241 17 18.942928 18.955540 3.303276e+07 \n", + "38242 17 18.921795 18.934738 3.303276e+07 \n", + "38243 18 18.893675 18.913025 3.303277e+07 \n", + "38244 20 18.867090 18.901897 3.303278e+07 \n", + "38245 20 18.852695 18.881191 3.303279e+07 \n", + "\n", + " layer_flag height_segment_ssh_flag geometry \\\n", + "0 0 0 POINT (-166.17149 -66.1212) \n", + "1 0 0 POINT (-166.17152 -66.12129) \n", + "2 0 0 POINT (-166.17154 -66.12139) \n", + "3 0 0 POINT (-166.17156 -66.12148) \n", + "4 0 0 POINT (-166.17159 -66.12156) \n", + "... ... ... ... \n", + "38241 0 0 POINT (18.94094 -63.78423) \n", + "38242 0 0 POINT (18.94093 -63.78418) \n", + "38243 0 0 POINT (18.94092 -63.78413) \n", + "38244 0 0 POINT (18.9409 -63.78406) \n", + "38245 0 0 POINT (18.94089 -63.784) \n", + "\n", + " hist_mean_median_h_diff \n", + "0 0.001873 \n", + "1 -0.005192 \n", + "2 -0.019470 \n", + "3 -0.045704 \n", + "4 -0.033787 \n", + "... ... \n", + "38241 -0.012611 \n", + "38242 -0.012943 \n", + "38243 -0.019350 \n", + "38244 -0.034807 \n", + "38245 -0.028496 \n", + "\n", + "[38246 rows x 12 columns]" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "gdf" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ce92fd39", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "2df0a0fb", + "metadata": {}, + "source": [ + "### Optical imagery to label point clouds\n", + "\n", + "Let's use the Sentinel-2 satellite image we found to label each ATL07 photon. We'll\n", + "make a new column called `sea_ice_label` that can have either of these classifications:\n", + "\n", + "0. thick/snow-covered sea ice\n", + "1. thin/gray sea ice\n", + "2. open water\n", + "\n", + "These labels will be empirically determined based on the brightness (or surface\n", + "reflectance) of the Sentinel-2 pixel.\n", + "\n", + "```{note}\n", + "Sea ice can move very quickly in minutes, so while we've tried our best to find a\n", + "Sentinel-2 image that was captured at about the same time as the ICESat-2 track, it is\n", + "very likely that we will still be misclassifying some of the points below because the\n", + "image and point clouds are not perfectly aligned.\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "c175a822", + "metadata": { + "lines_to_next_cell": 2 + }, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + "\n", + "\n", + "
\n", + "
\n", + " \n", + "
\n", + "
" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Get first STAC item from collection\n", + "item = item_collection.items[0]\n", + "item" + ] + }, + { + "cell_type": "markdown", + "id": "ee3e932a", + "metadata": {}, + "source": [ + "Use [`stackstac.stack`](https://stackstac.readthedocs.io/en/v0.5.1/api/main/stackstac.stack.html)\n", + "to get the RGB bands from the Sentinel-2 image." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "7e6b8e62", + "metadata": { + "lines_to_next_cell": 2 + }, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'stackstac-bd9892c19b6c86ec2d80ee484fd3c349' (band: 3,\n",
+       "                                                                y: 1831, x: 1830)> Size: 40MB\n",
+       "dask.array<getitem, shape=(3, 1831, 1830), dtype=float32, chunksize=(1, 1024, 1024), chunktype=numpy.ndarray>\n",
+       "Coordinates: (12/54)\n",
+       "    time                                     datetime64[ns] 8B 2019-10-31T20:...\n",
+       "    id                                       <U23 92B 'S2B_2CND_20191031_0_L2A'\n",
+       "  * band                                     (band) <U5 60B 'red' 'green' 'blue'\n",
+       "  * x                                        (x) float64 15kB 5e+05 ... 6.097...\n",
+       "  * y                                        (y) float64 15kB 1.9e+06 ... 1.7...\n",
+       "    s2:sequence                              <U1 4B '0'\n",
+       "    ...                                       ...\n",
+       "    title                                    (band) <U20 240B 'Red (band 4) -...\n",
+       "    raster:bands                             object 8B {'nodata': 0, 'data_ty...\n",
+       "    common_name                              (band) <U5 60B 'red' 'green' 'blue'\n",
+       "    center_wavelength                        (band) float64 24B 0.665 0.56 0.49\n",
+       "    full_width_half_max                      (band) float64 24B 0.038 ... 0.098\n",
+       "    epsg                                     int64 8B 32702\n",
+       "Attributes:\n",
+       "    spec:        RasterSpec(epsg=32702, bounds=(499980, 1790160, 609780, 1900...\n",
+       "    crs:         epsg:32702\n",
+       "    transform:   | 60.00, 0.00, 499980.00|\\n| 0.00,-60.00, 1900020.00|\\n| 0.0...\n",
+       "    resolution:  60
" + ], + "text/plain": [ + " Size: 40MB\n", + "dask.array\n", + "Coordinates: (12/54)\n", + " time datetime64[ns] 8B 2019-10-31T20:...\n", + " id " + ] + }, + "metadata": { + "image/png": { + "width": 500 + } + }, + "output_type": "display_data" + } + ], + "source": [ + "fig = pygmt.Figure()\n", + "\n", + "# Plot Sentinel-2 RGB image\n", + "# Convert from 14-bit to 8-bit color scale for PyGMT\n", + "fig.grdimage(grid=((da_image / 2**14) * 2**8).astype(np.uint8), frame=True)\n", + "\n", + "# Plot ATL07 points\n", + "# Sea ice points in blue\n", + "df_ice = gdf[gdf.height_segment_ssh_flag == 0].get_coordinates()\n", + "fig.plot(x=df_ice.x, y=df_ice.y, style=\"c0.2c\", fill=\"blue\", label=\"Sea ice\")\n", + "# Sea surface points in orange\n", + "df_water = gdf[gdf.height_segment_ssh_flag == 1].get_coordinates()\n", + "fig.plot(x=df_water.x, y=df_water.y, style=\"c0.2c\", fill=\"orange\", label=\"Sea surface\")\n", + "fig.legend()\n", + "fig.show()" + ] + }, + { + "cell_type": "markdown", + "id": "294f2407", + "metadata": { + "lines_to_next_cell": 2 + }, + "source": [ + "Looking good! Notice how the orange (sea surface) coincide with the cracks in the sea\n", + "ice." + ] + }, + { + "cell_type": "markdown", + "id": "7655a621", + "metadata": {}, + "source": [ + "Next, we'll want to do better than a binary 0/1 or ice/water classification, and pick\n", + "up different shades of gray (white thick ice, gray thin ice, black water). Let's\n", + "use the X and Y coordinates from the point cloud to pick up surface reflectance values\n", + "from the Sentinel-2 image.\n", + "\n", + "PyGMT's [`grdtrack`](https://www.pygmt.org/v0.12.0/api/generated/pygmt.grdtrack.html)\n", + "function can be used to do this X/Y sampling from a 1-band grid." + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "e6779b6c", + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "grdtrack [WARNING]: Some input points were outside the grid domain(s).\n" + ] + } + ], + "source": [ + "df_red = pygmt.grdtrack(\n", + " grid=da_image.sel(band=\"red\").compute(), # Choose only the Red band\n", + " points=gdf.get_coordinates(), # x/y coordinates from ATL07\n", + " newcolname=\"red_band_value\",\n", + " interpolation=\"n\", # nearest neighbour\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "f457a8f0", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk0AAAHFCAYAAADv8c1wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAACwS0lEQVR4nOzdeVhUZfsH8O+wC8LIjggCKi64JLkgiluuCfIz39TSsMWl0jQNFyrzNTPFLSuXyjRNK+W10lBIQVMUBReEzA1BdgVZRBBQ1vP7g2ZkmJlzzuwzcH+uy6s455mZZ2CW+zzLfQsYhmFACCGEEEJYGem6A4QQQgghhoCCJkIIIYQQHihoIoQQQgjhgYImQgghhBAeKGgihBBCCOGBgiZCCCGEEB4oaCKEEEII4YGCJkIIIYQQHihoIoQQQgjhgYImolcuXryIl156CR07doS5uTmcnZ3h7++P0NBQjT5uVVUVVq1ahTNnzkid27t3LwQCAbKysjTaB09PT7zxxhuc7fLz87FixQr4+/vDwcEBNjY26NevH3bu3In6+nqN9pHLqlWrIBAIONuNGDECvXr10kKP+NPW31lfHpcvvq9LXRsxYgRGjBgh/jkrKwsCgQB79+7Vel+io6OxatUqmecM5fdJZDPRdQcIEYmKikJwcDBGjBiBDRs2oH379sjPz8eVK1dw8OBBbN68WWOPXVVVhU8//RQAJD54ASAwMBAJCQlo3769xh5fEUlJSdi3bx9mzpyJTz75BKampvjzzz/x7rvvIjExET/88IOuu0hakMOHD8PGxkbX3VBY+/btkZCQgM6dO2v9saOjo7F9+3aZgZOh/j5JIwqaiN7YsGEDvLy8cOLECZiYPHtpvvLKK9iwYYPO+uXo6AhHR0edPX5zQ4YMwd27d2Fqaio+NmbMGNTU1GD79u349NNP4e7urtB91tbWQiAQSPzeCQEAX19fXXdBKebm5hg0aBBnu6qqKlhaWmqhR40M9fdJGtH0HNEbJSUlcHBwkPnFbWQk/VKNiIiAv78/rKys0LZtW4wbNw7JyckSbd544w20bdsW6enpmDBhAtq2bQt3d3eEhoaiuroaQOMwvigo+vTTTyEQCCAQCMRD6LKmT0TTS5cvX8bQoUNhaWmJTp06ITw8HA0NDRJ9KC8vx5IlS+Dl5QUzMzN06NABixYtQmVlpVK/J1tbW4mASWTgwIEAgLy8PNbbnzlzBgKBAPv370doaCg6dOgAc3NzpKenAwBOnjyJUaNGwcbGBpaWlhgyZAhOnToldT9RUVHo27cvzM3N4eXlhU2bNin8XM6dO4dBgwahTZs26NChAz755BOpKcZPP/0Ufn5+sLOzg42NDZ5//nns3r0bzWuNe3p6IigoCMePH8fzzz+PNm3aoHv37jJH3hITEzFkyBBYWFjA1dUVH374IWprazn7++WXX0IgEIh/V00tX74cZmZmKC4uBgDExsbi//7v/+Dm5gYLCwt06dIFb7/9tvg8G3lTOM2noAD+r69Dhw7Bz88PQqFQ/Hp96623FO6L6PVz4MABfPzxx3B1dYWNjQ1Gjx6N1NRU1vs6d+6c+LbN7du3DwKBAJcvX2a9D4ZhsGHDBnh4eMDCwgLPP/88/vzzT6l2sqbnRNPHV69excsvvwxbW1vxSBTDMNixYwf69u2LNm3awNbWFi+//DIyMjKk7vv48eMYNWqU+HfZo0cPrFu3DkDjZ8727dsBQPxZ0vTzQ9bfNicnB6+99hqcnJxgbm6OHj16YPPmzRKfJaLns2nTJnzxxRfw8vJC27Zt4e/vj8TERNbfGVEjhhA9MXv2bAYAs2DBAiYxMZGpqamR2/bzzz9nBAIB89ZbbzHHjh1jfv/9d8bf35+xsrJibty4IW73+uuvM2ZmZkyPHj2YTZs2MSdPnmRWrlzJCAQC5tNPP2UYhmGePn3KHD9+nAHAzJo1i0lISGASEhKY9PR0hmEYZs+ePQwAJjMzU3y/w4cPZ+zt7Rlvb2/m22+/ZWJjY5l58+YxAJgff/xR3K6yspLp27cv4+DgwHzxxRfMyZMnma+++ooRCoXMCy+8wDQ0NIjbenh4MK+//rrSv7/XX3+dMTExYYqLi1nbnT59mgHAdOjQgXn55ZeZyMhI5tixY0xJSQmzf/9+RiAQMJMmTWJ+//135ujRo0xQUBBjbGzMnDx5UnwfJ0+eZIyNjZmAgADm999/Zw4dOsQMGDCA6dixI8PnY0X0+3N1dWW+/vpr5sSJE8zChQsZAMz8+fMl2r7xxhvM7t27mdjYWCY2Npb57LPPmDZt2oj/fiIeHh6Mm5sb4+Pjw+zbt485ceIEM2XKFAYAExcXJ25348YNxtLSkvHx8WEOHDjA/PHHH8y4cePEfW/6d26uqKiIMTMzYz7++GOJ43V1dYyrqyszefJk8bFvvvmGWbduHRMZGcnExcUxP/74I/Pcc88x3bp1k3hty3p9yXstDB8+nBk+fLj4Z76vrwsXLjACgYB55ZVXmOjoaOavv/5i9uzZw4SEhMh9rvL6Inr9eHp6MjNmzGCioqKYAwcOMB07dmS8vb2Zuro61vvz9fVlhgwZInV8wIABzIABAzj789///lf8Xv3zzz+ZnTt3Mh06dGBcXFwkfjeZmZkMAGbPnj1St/Xw8GCWL1/OxMbGMkeOHGEYhmHmzJnDmJqaMqGhoczx48eZX375henevTvj7OzMFBQUiO9j165djEAgYEaMGMH88ssvzMmTJ5kdO3Yw8+bNYxiGYdLT05mXX36ZASD+LElISGCePn0q8/dZWFjIdOjQgXF0dGS+/fZb5vjx48x7773HAGDeffddqefj6enJjB8/njly5Ahz5MgRpnfv3oytrS3z6NEjzt8dUR0FTURvFBcXMwEBAQwABgBjamrKDB48mFm3bh3z+PFjcbucnBzGxMSEWbBggcTtHz9+zLi4uDBTp04VH3v99dcZAMz//vc/ibYTJkxgunXrJv65qKiIAcD897//leqXvKAJAHPx4kWJtj4+Psy4cePEP69bt44xMjJiLl++LNHu119/ZQAw0dHR4mOqBE0nTpxgjIyMmMWLF3O2FX3pDRs2TOJ4ZWUlY2dnx0ycOFHieH19PfPcc88xAwcOFB/z8/NjXF1dmSdPnoiPlZeXM3Z2dryDJgDMH3/8IXF8zpw5jJGREZOdnS3zdvX19UxtbS2zevVqxt7eXirotLCwkLjtkydPGDs7O+btt98WH5s2bRrTpk0biS/Curo6pnv37pxBE8MwzOTJkxk3Nzemvr5efCw6OpoBwBw9elTmbRoaGpja2lomOztb6nmrEjTxfX1t2rSJAaDUF6u8oGnChAkS7f73v/+JAwU2ouebnJwsPnbp0iWpCw5ZSktLGQsLC+all16SOH7+/HkGAO+gaeXKlRK3T0hIYAAwmzdvljiem5vLtGnThlm2bBnDMI2fMTY2NkxAQIDEa6+5+fPny30fNP99hoWFyfwseffddxmBQMCkpqZKPJ/evXtLBKai392BAwfk9oeoD03PEb1hb2+Pc+fO4fLlywgPD8f//d//4c6dO/jwww/Ru3dv8bTGiRMnUFdXh5kzZ6Kurk78z8LCAsOHD5faAScQCDBx4kSJY3369EF2drZK/XVxcRFPicm732PHjqFXr17o27evRF/HjRsHgUAgc7eeSH19vcRtmk/7iVy9ehVTp07FoEGDxFMEfPznP/+R+PnChQt4+PAhXn/9danHHT9+PC5fvozKykpUVlbi8uXLmDx5MiwsLMS3t7a2lvo9s7G2tkZwcLDEsenTp6OhoQFnz54VH/vrr78wevRoCIVCGBsbw9TUFCtXrkRJSQkKCwslbt+3b1907NhR/LOFhQW6du0q8Tc5ffo0Ro0aBWdnZ/ExY2NjTJs2jVe/33zzTeTl5eHkyZPiY3v27IGLiwtefPFF8bHCwkK88847cHd3h4mJCUxNTeHh4QEAuHXrFq/H4sL39TVgwAAAwNSpU/G///0P9+7dU/mxm//t+vTpAwCc76tXX30VTk5O4iksANi6dSscHR3Ff4OGhgaJ5yOask1ISMDTp08xY8YMifscPHiw+HfLR/PX/rFjxyAQCPDaa69JPK6Liwuee+458e/xwoULKC8vx7x583jtEuXjr7/+go+Pj9RnyRtvvAGGYfDXX39JHA8MDISxsbH4Z76/d6IeFDQRvdO/f38sX74chw4dwv3797F48WJkZWWJF4M/ePAAQOMXgampqcS/iIgIqTUjlpaWEl/uQOMi0adPn6rUT3t7e6lj5ubmePLkifjnBw8e4Nq1a1L9tLa2BsMwrOtbOnfuLHGb1atXS7VJTk7GmDFj4O3tjejoaJibm/Puf/PdgKLf68svvyzV3/Xr14NhGDx8+BClpaVoaGiAi4uL1H3KOiZP06Cl+e1LSkoAAJcuXcLYsWMBAN9//z3Onz+Py5cv4+OPPwYAid81wO9vUlJSolLfX3zxRbRv3x579uwBAJSWliIyMhIzZ84Uf5k1NDRg7Nix+P3337Fs2TKcOnUKly5dEq89ad5vZfF9fQ0bNgxHjhwRX2y4ubmhV69eMtcW8dX8dy167XE9N3Nzc7z99tv45Zdf8OjRIxQVFeF///sfZs+eLb6P1atXSzwf0boj0etC1deerNc+wzBwdnaW+l0mJiaKf49FRUUAADc3N96PxaWkpETmzlxXV1fx+aaU/b0T9aCtMkSvmZqa4r///S+2bNmC69evAwAcHBwAAL/++qtCV5e64ODggDZt2shNAyB6LrIcPXpUvFgdePYhKpKcnIzRo0fDw8MDMTExEAqFCvWt+ZWyqC9bt26Vu+vI2dlZvNOuoKBA6rysY/KIgjRZtxd9MRw8eBCmpqY4duyYROB75MgR3o/TnL29vUp9NzY2RkhICL7++ms8evQIv/zyC6qrq/Hmm2+K21y/fh1///039u7di9dff118XNYCclksLCwk/vYixcXFEq8ZRV5f//d//4f/+7//Q3V1NRITE7Fu3TpMnz4dnp6e8Pf359UvdXn33XcRHh6OH374AU+fPkVdXR3eeecd8fm5c+ciKChI/LMoMBC9LuT9/Tw9PXk9vqzXvkAgwLlz52ReeIiOiTaMcG22UIS9vT3y8/Oljt+/f1/cN6I/KGgieiM/P1/mFZdoKkMUNIwbNw4mJia4e/eu1DC7sjR1tRYUFIS1a9fC3t4eXl5eCt22d+/ecs+lpKRg9OjRcHNzQ2xsLGxtbVXtKoYMGYJ27drh5s2beO+99+S2MzMzw8CBA/H7779j48aN4mDm8ePHOHr0KO/He/z4MSIjIyWmeX755RcYGRlh2LBhACBOg9B0OuLJkyfYv3+/ok9PbOTIkYiMjMSDBw/Eo1319fWIiIjgfR9vvvkmNmzYgAMHDmDv3r3w9/dH9+7dxedFX8rNv4C/++47Xvfv6emJa9euSRy7c+cOUlNTJb5ElXl9mZubY/jw4WjXrh1OnDiB5ORkrQdN7du3x5QpU7Bjxw7U1NRg4sSJEtOqrq6uUhcJADBo0CBYWFjg559/lnjvX7hwAdnZ2byDpuaCgoIQHh6Oe/fuYerUqXLbDR48GEKhEN9++y1eeeUVuVN0TT9P2rRpw/rYo0aNwrp163D16lU8//zz4uOi3YQjR45U4hkRTaGgieiNcePGwc3NDRMnTkT37t3R0NCAlJQUbN68GW3btsX7778PoPELZfXq1fj444+RkZGB8ePHw9bWFg8ePMClS5dgZWUlTlTJl7W1NTw8PPDHH39g1KhRsLOzg4ODg9IfwiKLFi3Cb7/9hmHDhmHx4sXo06cPGhoakJOTg5iYGISGhsLPz0+h+0xNTcXo0aMBAJ9//jnS0tKQlpYmPt+5c2el8kq1bdsWW7duxeuvv46HDx/i5ZdfhpOTE4qKivD333+jqKgI33zzDQDgs88+w/jx4zFmzBiEhoaivr4e69evh5WVFR4+fMjr8ezt7fHuu+8iJycHXbt2RXR0NL7//nu8++674i/QwMBAfPHFF5g+fTrmzp2LkpISbNq0SaFpyOZWrFiByMhIvPDCC1i5ciUsLS2xfft2hVJAdO/eHf7+/li3bh1yc3Oxc+dOqfOdO3dGWFgYGIaBnZ0djh49itjYWF73HxISgtdeew3z5s3Df/7zH2RnZ2PDhg1Sf1e+r6+VK1ciLy8Po0aNgpubGx49eoSvvvoKpqamGD58OO/nrU7vv/+++LUvmurkYmtriyVLlmDNmjWYPXs2pkyZgtzcXKxatUqh6bnmhgwZgrlz5+LNN9/ElStXMGzYMFhZWSE/Px/x8fHo3bs33n33XbRt2xabN2/G7NmzMXr0aMyZMwfOzs5IT0/H33//jW3btgF4dsGzfv16vPjiizA2NkafPn1gZmYm9diLFy/Gvn37EBgYiNWrV8PDwwNRUVHYsWMH3n33XXTt2lXp50U0QJer0AlpKiIigpk+fTrj7e3NtG3bljE1NWU6duzIhISEMDdv3pRqf+TIEWbkyJGMjY0NY25uznh4eDAvv/yyxNb4119/nbGyspK6rWgXTVMnT55kfH19GXNzcwaAeIeLvN1zPXv2lLrf119/nfHw8JA4VlFRwaxYsYLp1q0bY2ZmxgiFQqZ3797M4sWLJXZw8d09J+qPvH9NdwvJItr9dOjQIZnn4+LimMDAQMbOzo4xNTVlOnTowAQGBkq1j4yMZPr06cOYmZkxHTt2ZMLDw2X+XmUR/f7OnDnD9O/fnzE3N2fat2/PfPTRR0xtba1E2x9++IHp1q0bY25uznTq1IlZt24ds3v3bpk7zgIDA2U+VtNdVQzTuNtq0KBBjLm5OePi4sIsXbqU2blzJ6/dcyKi9m3atGHKysqkzt+8eZMZM2YMY21tzdja2jJTpkxhcnJypHZpynp9NTQ0MBs2bGA6derEWFhYMP3792f++usvmc+Fz+vr2LFjzIsvvsh06NCBMTMzY5ycnJgJEyYw586d43ye8nbPNX89yNqtxsXT05Pp0aMH7/YM0/i7WbduHePu7s6YmZkxffr0YY4ePSr1u2HbPVdUVCTzvn/44QfGz8+PsbKyYtq0acN07tyZmTlzJnPlyhWJdtHR0czw4cMZKysrcfqK9evXi89XV1czs2fPZhwdHRmBQCDx95X1Ps/OzmamT5/O2NvbM6ampky3bt2YjRs3SuzQFD2fjRs3SvW7+WuKaI6AYZpliCOEEEI07Nq1a3juueewfft2zJs3T9fdIYQXCpoIIYRozd27d5GdnY2PPvoIOTk5SE9P12oZE0JUQSkHCCGEaM1nn32GMWPGoKKiAocOHaKAiRgUGmkihBBCCOGBRpoIIYQQQnigoIkQQgghhAcKmgghhBBCeKDklmrU0NCA+/fvw9raWm3FHAkhhBCiWQzD4PHjx3B1dYWRkfzxJAqa1Oj+/ftwd3fXdTcIIYQQooTc3FzWgswUNKmRtbU1gMZfuo2NjY57QwghhBA+ysvL4e7uLv4el4eCJjUSTcnZ2NhQ0EQIIYQYGK6lNTpdCH727FlMnDgRrq6uEAgEOHLkiNy2b7/9NgQCAb788kuJ49XV1ViwYAEcHBxgZWWF4OBg5OXlSbQpLS1FSEgIhEIhhEIhQkJC8OjRI4k2OTk5mDhxIqysrODg4ICFCxeipqZGTc+UEEIIIYZOp0FTZWUlnnvuOXFlaHmOHDmCixcvwtXVVercokWLcPjwYRw8eBDx8fGoqKhAUFAQ6uvrxW2mT5+OlJQUHD9+HMePH0dKSgpCQkLE5+vr6xEYGIjKykrEx8fj4MGD+O233xAaGqq+J0sIIYQQw6bDYsESADCHDx+WOp6Xl8d06NCBuX79OuPh4cFs2bJFfO7Ro0eMqakpc/DgQfGxe/fuMUZGRszx48cZhmmsNA6ASUxMFLdJSEhgADC3b99mGKaxYrWRkRFz7949cZsDBw4w5ubmMquXy1NWVsYAUOg2hBBCCNEtvt/fep2nqaGhASEhIVi6dCl69uwpdT4pKQm1tbUYO3as+Jirqyt69eqFCxcuAAASEhIgFArh5+cnbjNo0CAIhUKJNr169ZIYyRo3bhyqq6uRlJQkt3/V1dUoLy+X+EcIIYSQlkmvg6b169fDxMQECxculHm+oKAAZmZmsLW1lTju7OyMgoICcRsnJyep2zo5OUm0cXZ2ljhva2sLMzMzcRtZ1q1bJ14nJRQKKd0AIYQQ0oLpbdCUlJSEr776Cnv37lU4USTDMBK3kXV7Zdo09+GHH6KsrEz8Lzc3V6F+EkIIIcRw6G3QdO7cORQWFqJjx44wMTGBiYkJsrOzERoaCk9PTwCAi4sLampqUFpaKnHbwsJC8ciRi4sLHjx4IHX/RUVFEm2ajyiVlpaitrZWagSqKXNzc3F6AUozQAghhLRsehs0hYSE4Nq1a0hJSRH/c3V1xdKlS3HixAkAQL9+/WBqaorY2Fjx7fLz83H9+nUMHjwYAODv74+ysjJcunRJ3ObixYsoKyuTaHP9+nXk5+eL28TExMDc3Bz9+vXTxtMlhBBCiJ7TaXLLiooKpKeni3/OzMxESkoK7Ozs0LFjR9jb20u0NzU1hYuLC7p16wYAEAqFmDVrFkJDQ2Fvbw87OzssWbIEvXv3xujRowEAPXr0wPjx4zFnzhx89913AIC5c+ciKChIfD9jx46Fj48PQkJCsHHjRjx8+BBLlizBnDlzaPSIEEIIIQB0PNJ05coV+Pr6wtfXFwDwwQcfwNfXFytXruR9H1u2bMGkSZMwdepUDBkyBJaWljh69CiMjY3FbX7++Wf07t0bY8eOxdixY9GnTx/s379ffN7Y2BhRUVGwsLDAkCFDMHXqVEyaNAmbNm1S35MlhBBCiEETMAzD6LoTLUV5eTmEQiHKyspohIoYpIyiCmQ/rIKnvRW8HKx03R1CCNEKvt/fVHuOEIJHVTVYeCAFZ9OKxMeGeTti66u+EFqa6rBnhBCiP/R2ITghRHsWHkjB+fRiiWPn04ux4ECyjnpECCH6h4ImQlq5jKIKnE0rQn2zmfp6hsHZtCJkFlfqqGeEEKJfKGgipJXLfljFej6rhIImQggBKGgipNXzsLNkPe9pTwvCCSEEoKCJkFavk2NbDPN2hHGzkkHGAgGGeTvSLjpCCPkXBU2EEGx91RdDujhIHBvSxQFbX/XVUY8IIUT/UMoBQgiElqbYN2sgMosrkVVSSXmaCCFEBgqaCCFiXg4ULBFCiDw0PUcIIYQQwgMFTYQQQgghPFDQRAghhBDCAwVNhBBCCCE8UNBECCGEEMIDBU2EEEIIITxQ0EQIIYQQwgMFTYQQQgghPFDQRAghhBDCAwVNhBBCCCE8UNBECCGEEMIDBU2EEEIIITxQ0EQIIYQQwgMFTYQQQgghPFDQRAghhBDCAwVNhBBCCCE8UNBECCGEEMIDBU2EEEIIITxQ0EQIIYQQwgMFTYQQQgghPFDQRAghhBDCAwVNhBBCCCE8UNBECCGEEMIDBU2EEEIIITxQ0EQIIYQQwgMFTYQQQgghPFDQRAghhBDCAwVNhBBCCCE8UNBECCGEEMKDToOms2fPYuLEiXB1dYVAIMCRI0fE52pra7F8+XL07t0bVlZWcHV1xcyZM3H//n2J+6iursaCBQvg4OAAKysrBAcHIy8vT6JNaWkpQkJCIBQKIRQKERISgkePHkm0ycnJwcSJE2FlZQUHBwcsXLgQNTU1mnrqhBBCCDEwOg2aKisr8dxzz2Hbtm1S56qqqnD16lV88sknuHr1Kn7//XfcuXMHwcHBEu0WLVqEw4cP4+DBg4iPj0dFRQWCgoJQX18vbjN9+nSkpKTg+PHjOH78OFJSUhASEiI+X19fj8DAQFRWViI+Ph4HDx7Eb7/9htDQUM09eUIIIYQYFkZPAGAOHz7M2ubSpUsMACY7O5thGIZ59OgRY2pqyhw8eFDc5t69e4yRkRFz/PhxhmEY5ubNmwwAJjExUdwmISGBAcDcvn2bYRiGiY6OZoyMjJh79+6J2xw4cIAxNzdnysrKeD+HsrIyBoBCtyGEEEKIbvH9/jaoNU1lZWUQCARo164dACApKQm1tbUYO3asuI2rqyt69eqFCxcuAAASEhIgFArh5+cnbjNo0CAIhUKJNr169YKrq6u4zbhx41BdXY2kpCS5/amurkZ5ebnEP0IIIYS0TAYTND19+hRhYWGYPn06bGxsAAAFBQUwMzODra2tRFtnZ2cUFBSI2zg5OUndn5OTk0QbZ2dnifO2trYwMzMTt5Fl3bp14nVSQqEQ7u7uKj1HQgghhOgvgwiaamtr8corr6ChoQE7duzgbM8wDAQCgfjnpv+vSpvmPvzwQ5SVlYn/5ebmcvaNEEIIIYZJ74Om2tpaTJ06FZmZmYiNjRWPMgGAi4sLampqUFpaKnGbwsJC8ciRi4sLHjx4IHW/RUVFEm2ajyiVlpaitrZWagSqKXNzc9jY2Ej8I4QQQkjLpNdBkyhgSktLw8mTJ2Fvby9xvl+/fjA1NUVsbKz4WH5+Pq5fv47BgwcDAPz9/VFWVoZLly6J21y8eBFlZWUSba5fv478/Hxxm5iYGJibm6Nfv36afIqEEEIIMRAmunzwiooKpKeni3/OzMxESkoK7Ozs4OrqipdffhlXr17FsWPHUF9fLx4NsrOzg5mZGYRCIWbNmoXQ0FDY29vDzs4OS5YsQe/evTF69GgAQI8ePTB+/HjMmTMH3333HQBg7ty5CAoKQrdu3QAAY8eOhY+PD0JCQrBx40Y8fPgQS5YswZw5c2j0iBBCCCGNtLGVT57Tp08zAKT+vf7660xmZqbMcwCY06dPi+/jyZMnzHvvvcfY2dkxbdq0YYKCgpicnByJxykpKWFmzJjBWFtbM9bW1syMGTOY0tJSiTbZ2dlMYGAg06ZNG8bOzo557733mKdPnyr0fCjlACGEEGJ4+H5/CxiGYXQSrbVA5eXlEAqFKCsroxEqQgghxEDw/f7W6zVNhBBCCCH6goImQgghhBAeKGgihBBCCOGBgiZCCCGEEB4oaCKEEEII4UGneZqI9o3aeBrZpVXwsrdCbOgIXXeHEEIIMRg00tRKfPrHP/AMi8LdkirUNQBpRZXwDIvC58du6LprhBBCiEGgoKmV2JOQI/P49/FZ2u0IIYQQYqAoaGoFRm08zXp+zOYz2ukIITJsO5WGV3cmYMfpdO7GhBCiQ7SmqRXILq1iPZ9ZUqmlnhDyzIX0Ikzf9ayQdkLGQ2w4kYqIOYPg19me5ZaEEKIbNNLUCnjYWrKe97K30lJPCHmmacDU1LTvE7XcE0II4YeCplbg1NKRrOdpFx3Rtm2n0ljP01QdIUQfUdDUSswJ8FToOCGadP5uMev5c2lFWuoJIYTwR2uaDNTqyBv49Woe6usbMKGPKzZOeY61/cdBPfFxUE+M2XwGmSWVlKeJ6NSQzg5IyHgo9/xQb0ct9oYQQvgRMAzD6LoTLUV5eTmEQiHKyspgY2OjkceIuZ6PuT9dlXnu62l9EezbQSOPS4i6eYZFyT2XFR6oxZ4QQlo7vt/fND1nYOQFTACwMCJFex0hREURcwYpdJwQQnSNgiYDsjqSO3v3i1+e0XxHCFEDv872yAoPxLJx3eDfyQ7LxnVDVnggpRsghOgtWtNkQC5ksC+eBYBbBZRziRiWeSO7YN7ILrruBiGEcKKRJgMyuJMDr3ZBX53VcE8IIYSQ1odGmgzIyuCe+OFCFme7O4UVvO4vo6gCx67dx8PKGozq4Uw7lghp4SZsiUN6cSW6OrXFsfeH6bo7hBgcCpoMjJuNOfLKq1nbdHVqy3r+UVUN5u67gktZpeJjey9kQ9jGBMfeGwp3e/YM4qT1GLXxNLJLqyhFhYFbH30T35zNFP98Pf8xPMOisGBkZ4SO667DnhFiWGh6zsDEfzSasw3XFeTCAykSAZNI2ZM6BG+PV7pvpOX49I9/4BkWhbslVahrANKKKuEZFoXPj3FvRiD6p2nA1NTW03e13BNCDBsFTQaoV3v5I0kv9W3PetuMogqcZcm2XFpVS9mYCfYk5Mg8/n18FhZFJOPQlVwt94goa8KWONbztAaSEP4oaDIg2SWV6LnyOK7ny1+zdDgln/0+HlZxPs7VHOlRKNJ6jNp4mvX8keT7WPrrNXT5KBo375VpqVdEWenF7Dtq+a6BJIRQ0GRQJm0/j8qaes52bMVOPey41ys939FWoX6RliW7lDuwBoC6BgbB289LHJuwJQ5dP46m0Qs9sT76Jmrq2Ys+cK2BBIAZOxPQe9VxhOxKVFfXCDFIFDQZiLjUQpRW1fJqK296bdTG0xjLMVRva2lKu+haOQ9b/hsB6hoYHLqSi/XRN+EZFoWbDypQU8+IFxpvPnFbgz0lXOStZWqKbQ3kzrh0eIZF4XzGQzx+Wo9z6SXwDIvCD+cy1NlNQgwG7Z4zECl5j3i3bR70fPrHP3LXqDTXp4NQkW4ZLNoVJl9np7a4W8JvtAkAdp69i7RC2VNAW0/fxZGreaiqb8AADzt8G9JfXd0kHLjWMgHAgpGdWc+v/TNV5vHVUbfw1tBOSvWLEENGI00Goq9bO95tm2dX5hswAcD1/Ja9RoV2hXFLylVsTZu8gEkkt6waJRW1OH7jATzDonDwYrYq3SM8ca1lMjMWyE03MGFLHGtBZQA0VUdaJQqaDMTwbk6wtTTlbNe82CnXot7mBnjYyT3XEtarsO0KI436uWt2TVvY4esavX/SqIuDFet5WWuZmk6zclFk9JuQloKCJgMSOT8AVmbGMs8F9XKRWeyU76JeEVnTJ1zrVUIjUhCw/hSWHvpbocfSNq4AcszmM9rpiJ57uZ8b77ayX43c3tl/RclbEr6iFw9nPS9rLROfNVAiiox+E9JS0JomA+Jub4kbq8fjXFoR1kXfwoPyp3ihuzM2TnlO7m08bC15r08Jf6mXzONsifGaJsc7lJSHQ0l5+HpaXwT7duD1mNrEFUBmlkhPZ2QUVSD0fynILKnEIC97hdfkGOLaqbk/XeXd1sgIqG9Q/DEuZz9U/EZEYQtGdpaZwFLWWiY+a6Ca2j97EHcjQloYAcMw7PtRCW/l5eUQCoUoKyuDjY2NrrsjxrU2QeT0khHwajak77vqT5Q+VfxbMSs8UOHbaNqojadZA0hvx2eBzaOqGrz8zQWkF0kHUuEv9cIrfh6sjyVv8f2cAE98HNRTsY6rCZ8A7uPfruHny5pPXDm+pzMtCteioK/O4k5hBWvNua4fR3OmJxBZGdiDFoKTFoXv9zcFTWqkr0HTO/su4/jNQl5t3WwsELVoGIT/rp/iG3A1N6WfG+sImK6wPZ+mgd7M3ZdYM6dzBYV8H0cbFAngeq48zisXmKp8nNuK1800DVaJ7kzYEse5lmloF/tWOcJkiCPGXETPydrcBOamxqirb8AkXzesCPLRddd0gu/3N61pagW+nTmA1yJyAMgrf4oFB5IBKD5c31RCRrHSt9WkOQGenMe5Ss0A7Gty9G3tFN/F71tibmslYAIg8eVMOxiVN2rjaXT5KEotrymuNVBZ4YEqB0xzf7wMz7AoeIZFYeBnMSrdlza0xN22zZ9T6ZM6FJRXo7iyFrviM+EZFoVTNx/oupt6i4KmViJyfgCMBfzank0rgveH/HbQyOPfyQFA4yLxPqtOoP9nMayZyvkKjUhBr5V/os9/j2PNsZsK3/7joJ7ICg+Et6MVTIwaRzmywgMlRlz4lJphW5OjzNopTVEkgPvqL8WKt5or+OnR1oz9Bvq2gzGjqAIfRKRg0vZ4tbx21UlTX+by8jZx5XPi8lNC45dxzK1nI96FlbXwDIvC4oP819BpW0vcbcsnBc2sfbRRQx6anlMjfZ2eE1F2qk0ZbkJz5JVVyzwXMWeQ1C4/LoeTcrH40DWZ53bP7I9RPs4K91GejKIKvLCZfZSNbU1OjxXReFIn/20lbzpqwpY4pBdXsq47UVSXj6JQx7IkzcQISF8biLk/Xpb4QmPTq701jr0/DP0/i0FxJb8s9QBgagTUciyP04epukdVNZj5w0VcyyuXOqfMa1cTND39y2cNlCK4PnsMfQ1kc5p4L/MVl1qI5b9eQ0lVNZ7r0A6/zhsiPsf1nJqaHeDVqqbqaHqOSFA0X5Oq5AVMADDte8mkeKsjb6DPqhPo+cmfctMWyAuYAPVfFXVybAtTjneGu20biZ9HbTwNr3+nHdgCJgBSH7SaLEHCVRLFy75x4T/fhJZmTYYrGZ4jlyJcAROg3VE4eRYeSJEZMAHSr11d0Mb077H3h+HO5xPU8oU/98fLnG3UMVUXl1qIoetPoccn0TITb87YmYAen/yJoetPyS011ZQyI8aKvJfVObUKNBZ0777iT7y+5zIKHlejth64kvMInmFR2P5XGq/n1FR8OvfvqDWilAMtxORt8bj14DF6treRuLIQUTRfkyKm9HPDoaQ8hW6z43Q6ujhaSW1vF6UtGN/TCTMGeWKotyNCI1I472/iV2dxVA0f8KEHr+K3lHzOdt/HZ+HjoJ4KlagBZK+pYkvpMLCTvUq1AE8tHcl6lS8K4Pq52/IaaWr6ReBpa4ES8B9p4kMUxOkKn/VsncKi4PPvaJsuaHv6d8KWOKQVVcBF2Ab7ZvlJ7bDlwicgL2QZsZywJQ63H1SgAYCVqRGMjQR4WtcgHsXJLqnE+C/P4kmTqFxUI29lYA/UNTRIlIPJLX2KkN2XYGEExIaOhLu99IXFhC1xrCO0wLPX6urIG7iQUYyALo7YFS//vSzKvt78M0M0tarqztpJ28/jqZxOb4y5A1NjAedzaiqgC9UglUWnI01nz57FxIkT4erqCoFAgCNHjkicZxgGq1atgqurK9q0aYMRI0bgxg3JOfvq6mosWLAADg4OsLKyQnBwMPLyJL/AS0tLERISAqFQCKFQiJCQEDx69EiiTU5ODiZOnAgrKys4ODhg4cKFqKmp0cTTVquvYlPhGRaFq3lleFLbIHVlIWJpquCwAE9fT+uLG0pkBj6XVsSaD+j4jUKE7L4E39UxOH+X+4rnn/zHCvdBFj4Bk8iYzWcUCphOLxmBPeez4BkWBe8PGwMZrsX2IbsvoftHUchVoBaciOiqV56mAdzO1wcofP9ZpU8Vvg0XcxPdDn7zWc/WAOi0IDHf0UNlzf3xMvqticHIjX+JR01qG4Dc0icYuekMAtaeQhnP4uEAvwzzTlbSG1WajtqIvusraxtQXl0vEbyP2XhGImBqanXULbn18542AMHb4+U+Jpel47o1Fi++kIXbBRVyAyYRUSUFVdZJyRuV51PQXd7vQZ7WNDWnCJ1+QlVWVuK5557Dtm3bZJ7fsGEDvvjiC2zbtg2XL1+Gi4sLxowZg8ePn31BLlq0CIcPH8bBgwcRHx+PiooKBAUFob7+2S6g6dOnIyUlBcePH8fx48eRkpKCkJAQ8fn6+noEBgaisrIS8fHxOHjwIH777TeEhoZq7smryZZTsheoboy5I/FzebV6l65N6eeGrPBABPt24KxxJUs1n7kaAKVVtXjIc92MMsPcTReW91gRrdBt02TkcGIzctMZiGbuahnw/nCW9eHeVFxqIabvTMSLX8ZJLFhWJLtzRpFyi/7b25gpdTt5dBmMAICHHXtA0pysxJEi8qZfVJmWibiUg94d27G2EY0eKlr2qOli7ZKKWmSWPJHZrukOWz74BORutpboviIaozadQWZxJT794x/er19VLm1Lq2olpur4PqaPi/QoOZc7hRUKTa2ujryBQWtPYtSm0/jvkWviAK38aR0qaxtwKCkPnmFRiEy+p/aSNj4u0iV2SCO9WQguEAhw+PBhTJo0CUDjKJOrqysWLVqE5cuXA2gcVXJ2dsb69evx9ttvo6ysDI6Ojti/fz+mTZsGALh//z7c3d0RHR2NcePG4datW/Dx8UFiYiL8/PwAAImJifD398ft27fRrVs3/PnnnwgKCkJubi5cXV0BAAcPHsQbb7yBwsJC3ou6tbkQnE9Olf4dGxcBTt4Wj6t56ivEa24MpH7+bOEmn740192lLW4XKL87TxbRomY+2BaW66v9swZKTNVll1Qi8OtzqKhWPk2AaAHu6dRCvLmHe+1Jc2bGAsx/oQu2xKZxN26iV3trXOcYHdTm4uCh4aeQ+6hx5MwYgCK/0V7NpurkTdn6uFjhZoF0oM1nWuafvEd4accF1DWwf1x3bGeBipo6PKyqkzpnYy7AHwuGy51eU3SjiKxkuPIcvJittzUHF4/xxvujuvL6HDMRAEZGAt5JQJvq1d4atx88Zp0iMxYA38x4XuGA7Mc3B+B1Jd6/8lhbGOOfVePVdn+GwOAXgmdmZqKgoABjx44VHzM3N8fw4cNx4cIFAEBSUhJqa2sl2ri6uqJXr17iNgkJCRAKheKACQAGDRoEoVAo0aZXr17igAkAxo0bh+rqaiQlJWn0eSpKkeHjG/mNi1lvPVDP1JVIdT0ktjdbWii+NE7dARMANDRAYsh6xs4EdF8Rjf6fxeDQFcks14YWMAHA1RzJtSGTtp9XKWACnl3ZKjrCItKhnYVSNci4AiYAWikM/dGvKfAMixIHTIBiARPQOIIQcSkHiyKScehKrtzpF1kBE9A4LeMVFoXOYVHo/GEUBq09icxmo7d8AiYAyHn0VGbABDSONo/cdAZdw6Kkptf4LNZuLkuBtVOv+Hno5Q45AHi+Y+P0IZ8R8zoGSgVMADCoswPn1Go9o1gZI5Fj1/J55+LjQxt1BflsANJHersQvKCgAADg7Cy5ldzZ2RnZ2dniNmZmZrC1tZVqI7p9QUEBnJycpO7fyclJok3zx7G1tYWZmZm4jSzV1dWorn62S6y8XPaOG3VSZMqlZ/vGaLmHs7VaR5pEvo/P0qt8JQ14tpC8qad1tVj66zUs/fUa/DvZwlWoXICga6IPd4DfGgY+RIuGOzm2xTBvR85F0NK3f4Lh3Zxga2mqlv40dadQ/YF1c79cuafyfdTUM1j++z8AgCPJ95W6Dwb/BmsMUFBejZGbzsDV2hx/Lh6O36/m8gqYePcXwHOrY2BlBqwK7oP/Xc7B5exHCt+Pp4Jrp7S9g5cv0ehtFwcrlXLTceFa86SKQ0l5GOBhi3/ul+GpjKUPS8d2lVqywUaRJKZzf7yMC3eL0dbCBKFju2NKf3fW9n6fx+DBY8nPCn2vW9qU3o40iQgEkguYGYaROtZc8zay2ivTprl169aJF5cLhUK4u7O/WFSlaIZu0S46v052muiOQUrIKMVvyap/UeqC6MP9q9hUtQ3FN100vPVVX3S0a8PSWrYdp9MROT8AFlx5GhTU1Umz6yqGhp/S6P2r6v7jaiw4kIxv4jI0cv+VNcDSX68pFTA5WZvjSpZiRZc1uYNXFV3+3ZTBlRFd313OLsXT2gbMGOgGVxtzmBo3LtHICg9Erw5CmPDcC7QysAevdk3XwFXUNKCgvAZLf21ce3XznvRFuigha/OAqamFPHZK65reBk0uLi4AIDXSU1hYKB4VcnFxQU1NDUpLS1nbPHggnRK+qKhIok3zxyktLUVtba3UCFRTH374IcrKysT/cnM1W+hU0QXXkf8GB4qMThH9NNzbQfz/8hb/K6NpziihpSnOLnsBp5eMgIct/+DpXFoR3O0tcfuzF9XWL6BxOkOTmk7J6auzaUUofCw/55muFD6uxtJfr6HLR9EyvyCbO5yUq9B2d21qmlbNro3efiXy9vOlPMwe2glpnwdi87S+8AyLwut7LoMjfRwAoI2JgHch5hV/yK/IELz9vNQxvjuNRWV2PMPUl8NKnfT2FeLl5QUXFxfExsaKj9XU1CAuLg6DBw8GAPTr1w+mpqYSbfLz83H9+nVxG39/f5SVleHSpUviNhcvXkRZWZlEm+vXryM//9l285iYGJibm6Nfv35y+2hubg4bGxuJf5rURcH8KAsjUlSqH0f0h7mJEXqt/BOd1ZjVXV4dPi8HK8QtfwE+zvxGepouTo+Yo75irpRcT//VNTAyvyCb0/c1hJ5hUTh4MRsVNXqxL0plq6NuYcSGvzB84xmFbmdlwW9dFNcauLoGRmIdqbJTs/pY60+nQVNFRQVSUlKQkpICoHHxd0pKCnJyciAQCLBo0SKsXbsWhw8fxvXr1/HGG2/A0tIS06dPBwAIhULMmjULoaGhOHXqFJKTk/Haa6+hd+/eGD16NACgR48eGD9+PObMmYPExEQkJiZizpw5CAoKQrdu3QAAY8eOhY+PD0JCQpCcnIxTp05hyZIlmDNnjl6VQ1Fm+FiTc/REe0RD4Ooqp9vDpS3nji2+r7d5I7uI/9+vsz2ywgMxVA2lRii5nmFo/gXZHJ/ktPog7PB1hS9M9VnWQ9kpI9hM6stvPRGfhKXn7z4r2q7q1Oz38Vl6EzjpNGi6cuUKfH194evrCwD44IMP4Ovri5UrVwIAli1bhkWLFmHevHno378/7t27h5iYGFhbW4vvY8uWLZg0aRKmTp2KIUOGwNLSEkePHoWxsbG4zc8//4zevXtj7NixGDt2LPr06YP9+/eLzxsbGyMqKgoWFhYYMmQIpk6dikmTJmHTpk1a+k3wp2jhTL7z2KR1uVVQITMJanNcrzd5I0v3Hyn+gd2cJpPrTd4mP+8VUVzTL8jmLmaVaLEnqlF2Z1xLwfc9xydh6ZAm0+tcuwb50JdNR3qTp6kl0GaepqCvzvLatq1MiRPSusjaCh6XWohPj95A7sMqmTXjBnm2w8F3pMv1iHiFRUGVDxZ1F2Furscnf8rNIk0Ut/HlPnJ3TY3edBrpxfq5CLw5EyPo7dorTVP0PceW18vESID0tRPEPytSKJhNW1MjXFfz2kkRg8/TRNgde38YrC2MOdttnPKcwqNTpHV5ecezNSnZJZXo8Ulj0c+MYtkBEwAkZj3CmmPyF4JamSn30WIM4NzSkRoNmIDGNBxEffbGS+/wC41IQcD6UwYTMAGAp13LmZ5TRFZ4oMLvufCXesk9Fzlf8oJKXTsnK2ob4BkWhT6rFKveoE4UNBmwPq5C1vPdnBs/AELHdW/cdtreGmbGAti10dv0XEQHRElQgcaEmXxHYHbFN245PnVTenfqJ0oUHrUwEWBQZ3vUaWHw+/f3AjT+GK3JjSbJag8n5cIzLAq/Jd9DngbqE2rSySUjdN0FnWC7AJLnFT8PfPRiN5nnEjMkp2TVMT3XVPlTBp5hUVj+a4pa75cPCpoM2M9z/VnPn1g8AttOpWH8l3GY8X0Clk/ogTufT8DV/47T2+y8RPtESVCVTZg5a98VqWPTBnZU+H6e1jE4f7cEIzedwbgv4hQqCquMpWO7avT+W5t39je+DvR9p5w80/o3LoJW5w5QQ6HsTlV5RYBXR92S+PnU0pFK3T+XCDUkp1UUBU0G7HCS/B0rr/m5wTMsCpti7+B2QQXO332IkN2X0GtlNHL/nVumaTsCPEuCqkrRT1lXqkrO0AEAUgsrFCoKq4z5L3jzTqtAuF3OfmgwO+VkEX0B/53HvTOspVFmp+qMnQms50N2JUr8bK3e2t5ig9ee1Mwdy0FBkwFju6L76aLsxd8VNQyCtzfuHBJN27nbWmikf0T/NR1tUaXelKwr1V4c08dczqYVSdVgUzdDzwKtTwZ42BnUTjlZlh76W+7oSUumzE7Va/fZk5o2vwhbEdRb4cfgI79cu8lfKWgyUKpc0ZVW1eLcvzXGfkrIRK6S6w46CDV06UA0zgiNiz/nv+AtPiaqIacMWVeq6lg3pEhRWGXRiKt6fBvSH36equfn0qU/UgyzxJIqds/sr9TtvDlyWjW/CEvI1ExA3d7GXCP3Kw8FTQZK1Su6qzmNQ9BsqfBF5gR4old7azRP+XSvrEalPhDdeSvAS+bxyPkBsDTj3pXZnLwrVVXXDSlaFFYZohFXoppOYVEwMtJsYjgnK1P4a7CWZmvL07RgZGeld6pmciTPbF70199LMwH1hY9Ga+R+5aGgyUCpekX3fEdbzlT4liYCZIUH4uOgnjj2/jDWvDtWppRF05DIC3Lc7S1xc/V47J81EF0crWBtYYzxPdk/VNmuVOe/4I2s8ED079gObRQs6DvM2xFeWsrQzLU+g3BrADSaE87ESIBLn4yVSJpIlLdgZGeEjuuu1G25No3MGCids0uZzSH6iIImA7V5Wl+Vbv9V7B3OVPi1DQz6rYnBO/uvcNYOchFaYJh3yy97wWcqx8QI8Ha0goYvupXm59WOtewF0FhP7mToCPyzajy+DekvNydL+Eu9eF2p/jpviEK5kbo5tcXWV315t1cV1/qM5vhWgmdjLBCgrbkxTJp9Co/v6SwRaDpbKTdlqm5CC8Ve0FP6ucHN1gLjezrj5ec7SEz99mpvDQVjaHHun/dGeXO0VE7/jqqtwTMUvdpbIys8UOmACeDeNOIklFwn+1VsKmsyTFUoky5BFZQRXI20mREcACKT72GhjLVNX0/ri9D/paCW5S/bxtQIQ7s4IOZWoVr6YmIEpK8NRGZxJbJKKrEm8gZnBlgzYwHc2rVBhhoyxWqDpakRbn72Imd2W29HK8SGjgDQeEW24EAyyp/Wic87WJliZDcnHLqqmfUT+2cNxOrIG0gralwPJABgaSpAZbMXhImRAJHzh8CnA78vi69iU7HlVLrU8aVju0qsjZKHbxbu/h3biXf0acuMnQk4n/GQV9vwl3rhFT8PAOxZkbkM83bE1ld9IeSxjuyrU3ewJZa95I25sQAnPxiBsV/GqZztPCs8EC/vOI8b+eXo2d5G/PdQ5PlO6eeGjVOeY23Ddn9tzYzQ1sIEoWO7S2Ubv3i3BNO+T5RzS8V0cbDEySUjEbD+FO+8UlnhgfBZEY2qOsP4+mz6maQOk7fF42qe/AuN5u9hTQVMANDdpS2OL1J9QwdlBG8Fgn07ICs8UHxF19PFGg5tzfDxkX/Q1oI9gWXP9jbYwPGBpgivf9eeeDlYYWQ3J155OWrqGYMJmAAgsI8rAO6cI7GhIyQyazcNmACguLIWh67eg+Irh/hZHXkDsaEjkBUeiKzwQGSGB6JaRqVfvhXqRWQFTACwMeYOr9vzHWnSdsAEcOc8E8kKDxQHTACwZIzia7be8O+I00tGYN+sgbwCJoDfzsYQf0+421vi1mcvYv+sgeho2wZtTAVwUHCkSrQO7dd5Q3Drsxcl/h5zAjx5309Chvx6dFz3NyfAE9dXv4jEj8bILM8iKgztadeGd3+Axr+ft6OVeDQ4KzwQJ5c0vp/5LnlwsW78fd5cM0H8Hutsr3ryxqZ9Uzd1BkwAkFrEPjKb+fDZBg5N13nUdmFvCppagJd8XZFX+hQ3Ch6juKIGj6vrUfqkjvU2v84bgoUHUtTWB1lvypa2K6npVTPbhz3AL7O2jDhGLTKb7TiLuJSDugbZV8RcFepFuD74mpZikYfPbjp1THspi+uxZb2eFZkqamtqhKzwQKz6v94Kr9Xis7Ox6Tq1od6OOLv8Bdz6bAIsFFjYv3tmf9ZRw4+Deoq/3Ln4d3q29igutRBD159Cj0+iJfL3NL2/poHMxzwyymcUVSCLYzFyU6K/X2zoCKSvDZT6zOK75CHx47FSx1RN3ihaFyjq27mlI6X+3gIAm6f0wfiezrAwEcCc559VkUCXr0qOXf4DPJ4t1v+bZURKHTRZ2FsWqqfRAkzfdUmh9kvHdkVGUQXOpimXBbY5eW9K0Zz51tN31fI48jhYmeDKJ+MAAGM2n0FmSaXai25+3ewD9eOgnvg4qKf48bzsJafklMmsLYuJkQAu1uYoePwU7dqYoLiSPRgGno36iXBt9T1/t1husVWRWw/Yi0M3LcXCZunYrjJHptzbWeBc2ChkFFVgwpdxyH5Yhec72krtwNGkt4Z2wltDOyFkVyJS8h6hvq4BdQzQ1aktjr0/TO7tIuYM4jVVpGqh0cj5ARj/9VlUyhg2ZFuM7+dpj7xS9qlgRadvRG3Zpl02TnkO2SWVGP/lWYkLiHPpJfAMi8LKwB54a2gniftTRPZD/qPUk3xdea3h+XpaX5lLHkTYLgTnBHji+/gsXv0xAeAsNMeLvV1lfum721sieeVYnEsrwtWcUjzf0RZD/10z+p9+ku/VHafT8f3Zu3haV4/hXZ1wt7BC6jNJnVZH3uBsY23xLOBT5uJwoKcd/s4rRTXH9GdnB/WWZ+GD1jSpkbbXNAHAtlNp2BTLb3pEZGgXe7w1tBPe3MO+e06W8T2dlXpTBn11FncKK9DVqS2yS6vw+Kn8t5K1hTHreRE/z3aIeEf2VM7Cn5MQ+U8B6+1nDfbA7gvZcs8bCYD/PM+9LqM5PutP+HCyMsWlT6SvakN2JeJcuvxAqPn2+YhLOVj++z9y27NVqBdRdA0Dl+brZR5V1WDyjvPIkFHctemXq7qFRqQg9tYDmBkL8FZAJ8wb2UWh26+OvIELGcXoaGeFGBk1+EQm9nHB1un9VO0uAOBcWhG2xNxBUcVTjOvZnteVNteaEmVTLrCtqwz27QDf1TGsFxCqpHrIKKrAC5vjeLXl8xpvaumhv3E6tRDlT2rQ0AB0d7FmDZybknUhFfTVWaQ+eAwXYRvsm+WntV2hmjD+yzjcblJrUBY3WwvELx+F0IgU/JbMvXazV3trbJ3RD1kllfC0txL/frw/jkIty1dBG1Mj3FLxYkSE7/c3BU1qpIug6dWdCUjguYBVxNrCGH/MD+D9gdOUOvLZcC26HdrFHoM62ctdK9O/oxC/zuOe6mH7wLa1NEXyyrHY/leazMfp2K4Nji4cynvNSVNxqYV4XYmAtDm23/Xnx27IvKqdE+Apc2qjy0fRMqfoTIwESF87gVd/2L54VX1dzNx9iXXkU915lA4n5crNqB8xZxD8OrOvb4m5no+5P13l/Xinl4zQ6RelvOAGeBbgqGLpob+RkFEM/04O4osMPu+DoV3sVRpN5HrdAIq9xgm31ZE38MOFLNY2ok0AfBbXD/CwxaF3B8s8p+6LNTa0ELyVUCZnSV+3dujk2FbhFAHytp3zNffHy+i3JgbWbdgDkf2zB2H+C94SiytdbMywf9ZAZIUH8gqYgMbpDBsZC+JtLIwROb/xPprmEbIwEaCzgxVOLxmBs2EvKBUwAYpl1jaXs+qTax2ComtBIucPgUmzHAii3XN8yUtUqWoCSz5Txc3rWClrxs4E9F51nLUEEZ+pNr4Bk7FAoNV8U/I03TTSxtQIbUwEmNLPDVnhgSoHTABwNesh8h89xbkmpW/41DJUpd4hAGx91ReDWQJcRV/jhNvKYO71ZqLAmc/i+pHdneSe41oHqYtNIzTSpEa6GGniisRl2T9rIIZ6O6Ksqhb/tz0eWTJ2sC0d2xX/3CvD5eyHGOBhh29DlEu1DzSWauGTeRxonIrJfViJPQk5UufkjaJwOZdWhMNX7wFg8NLzbuK1AZqUW1LFufXbXACcXDIS7vaWMof0NeHQlVycv1uMIZ0dFJquaErWVnRVnE4t5JwqtrYwxj+rxiv9GDvj0hWqKbZsXDe5U3V8rrRFFEkroG2jNp5GdmmVSq+5T//4R+Z71dXaHCsm+mDeL+xFl9lGmhTpX2ZxJS5mlIABUFpZgzuFj1V6jRN2p24+wKx9V2Seaz5yqerUsLzZAL6pTvii6Tkd0EXQxDf3TXNmaNzx4W5vicziSsz98TJySqvQ21Wo1uj9UVUN+q6OlXt+aBd7pOQ9Ql+3duIPT01OA6lqwpY43HzQOJ/fqz33OodzaUX4+Pd/UFTxFMYQoKaBgZ2VKTZO6auV4M0Q8Fmbouo0jqJ5Yvw72eGAnDQEXGs6zI0F+HZmf4m1GfpEXqCjzEUJ2+91mLcj/rn3SOE1TersH9GsNcdu4n9XclFf34AJfVxlrv/s/nEU2Jao8vkcBdR/sdYcBU06YCgjTSKidT2a1O+zWJRUyq9RN76ns8QoliKJI7VpffRNfHM2U+a55uUIJmyJQ3pxJeeuK/KMJtc0KZK4UoRtpCnk+wScu8t9f+q+ElYXtkBn98z+vGuRcb1XAeCX2X6Y9eNlmRd28hb46/NFE+GPbd1gU2bGAtz5XPdrzmhNUyuhSiX50qpanFMi7cDqyBvovyYWg9bGYsdp2QkPgcYRBLaACQBim+04yi5l/xBunoNIGy6kF8kNmIBnKRXWR9+EZ1gUbj6oQE09g+v5j+EZFoXNJ25rq6sGa+urvugsJ/ePqrmbFC2RAoB1Fx2fgAngn/RTm7jKIcmbcpGF670KANX1DVLJNod2aUxMKStg4urfmM1nePeP6BafgAloTOlhSChoagFUWYh7NYe9/lxTMdfz4RkWhR8uZKG4ogYF5TXYcKKxptDFu9Jb4PnkUWleVNzDlj3vRvMcRNrAJw9W0Fdn5QZWms5T1RIILU0xrb+bRu67j6tiNcUi5sifBgxlyeEjC5+kn9rEJ9DhW8uL670KAJ7/vl+bJttkm2bVx4smojhF3idco/ETtsSh68fRCPrqrIq9Ug+Vgqb09HScOHECT540ZmWlmT7daLoDTNEU/L9eyRPvduHCtmNI1o4jDzvFE4/xKVGiTdtO8cu3dD2fPfmjvrzh9Zm8hdqro26pdL9cJVIEAMyNG6fkssIDWdMNXMxiTxTaHN+kn9rCJ9CJT+c3+sz1XlVm16A+XjQRxZ26LT9nWVNsyUL1deReqaCppKQEo0ePRteuXTFhwgTk5+cDAGbPno3Q0FC1dpDw9+u8Ibw+FJvKLX2CkZvOYOSG0yhjWbDJJwts86m6To5t4Wxtznm7d/ZLTglwlSjRpvN3uetn8XGnkD0ZXGs3Y2cC63lVUw6wTfExAKrrgQ0nUnGKJUElwL8+mUjP9tpZ28gXn3IfitTykveedLU2x9ZXfXnfj4i+XTQR5VRWc1cuWDq2K0LHdceojafR5aMoqalXfR25VypoWrx4MUxMTJCTkwNLy2df0tOmTcPx48fV1jmiuPvl/Kp0N5f5sAoLDsjfHnyBR/FNWeuj+nnact7ucrbkGhFV6lGpG988WO7tLFjPG9q8vbZxrTtSNZ/PW0M7ISs8EEO72MPaQn7RLq41PZlFigW/usgjw4Xr4kORWl5N36vGAsDFxhynl4zAhY9H806zMGNnArqviEb/z2Jw6EquXl00tRai/GXqyofGZ85pY8wdeIZF4W5JFeoagLSiSniGReHzYzcwYQv7blpdjtwrVXsuJiYGJ06cgJub5BoEb29vZGfLL0tBNK+Hs7XSu+nO/puYTtaQ+uBODpyp82VtoR/u7YhojnImTYs7NqUPV5XvjfLmVabGwpS9eibtomPXx1XIusOtr1s7tTxOQwODJzXsJXrWHLspM3DILqlU6L2latJPTfk4qCcGdXKQGSCy1bAD5JedUea92jx31tO6Wiz9tXHxcPSCALx/MFkructas+Z/A1l1AZXRx1Wo9PfQ9/FZMDMWsLbR5ci9UiNNlZWVEiNMIsXFxTA3556OIZqjym46AMiSs9CSTxZYWTuOpg3sKJWJujlVEmdqA9vCYJHsh1Vy5+fZ5u1JI651R6oW7t0Zlw7PsCicz3jIWcxZ3pqeSdv5Leru37EdssID9TLdgMgoH2dkhQdidoAXuru0xewAL2SFB8pNN3A4KReeYVH4Lfkeyp/WobiylnUTCB9syUaDt59HbOgIpK8NpIBJgzS1jtCvk+wLYb4aZJR8akqXI/dKBU3Dhg3Dvn37xD8LBAI0NDRg48aNGDmSe86caNaswR2Vvq0ny0JLtqtQtsAicv4QyIubVC3Nog1+ne3hYMk+ktTVqS1Cx3VHVnggerW3hpmxAL3aWyMrPJBXdXVtEZWyab6OTNd2xslPXaFqygGA/Qu6OVOB9Is1LrWQNUmjSGOZH/2bkpNnRZAPji8azjklp2rZmea41rDVNTA4dCVX4fsl/GlyHSFbihZe2K+zdTpyr9T03MaNGzFixAhcuXIFNTU1WLZsGW7cuIGHDx/i/Hn92mLbGn0S3BtpRU8463k1x7XbRXR1uubYTRxJuQcTI2CmvxdnZXifDkJkrAvEoSu5+CImFeVPazHU21HvR5iaurJyPGvSvaZvYn2cimteyub4jQfwDItC+Eu98Iqfhw571ogtqFFlmgDg/nJo7h8ZOyH5rKlyaGum0OMYCj7bx3ecTuf8HGiKT+6s83eLqQyKBmliHaG6ipV72VthfC8XmYu+dT1yr9RIk4+PD65du4aBAwdizJgxqKysxOTJk5GcnIzOnWkqQh9sfdVXqiCvuYkRerpay2zvZWfJe7fLiiAfXFkxBokfjVHog3JKf3ckfDQaN1a/aFABk4ghT7/Jq/0Xdvi6lnsiLfhr9kWdqi5OvZKlWDZwQDqJIp81VZP6ql74Vh/xSbOgaJJcPrmzlClGTvjj+hsoso4wu6QSnT+MUkvABDSukdPXkXsqo6JGuiijwiWzuBJZJZVSdbBCdiUiOacUHe2ssP21fnpZI0tfBX11FncKKwymTMrcHy8j5lah3PPNS9loG1ddOFWL9Spadw4ATIyA9LWSJTt8V8coXEetJQiNSMFvyfc42zUvJ8SF7e9iYiRA+lrdl9Zo6dRRsoarvqiidFVjkO/3t1LTc2fPsl8ZDhum/18krYWXg+yioaourG3NDCFQaioplz3re/OUD9oUcUm6MGtzRipc1imawVtEVhLFyPkBGP/1WVRWS+++49p1Zsg2T+vLK2jaevquQkHTysAechccR843nHVhhkze34DvOsLm0/6KEBXqHbP5jEHtklRqpMnISHpWT9Bk8WR9PfuW3pZKH0eaCNHnkaZFEck4knyftY0qI00B608hr1Tx3GVsV9nn0oqwJeYOiiqeYlzP9grlNTJUkcn3sJBHAMq3Yn1TIbsScTnrIdqam2D5iz1oHZMOhOxKREreI/R1a6fQBbUyo7gi+jYyq9GRptJSySvX2tpaJCcn45NPPsHnn3+uzF0SAzVq42lkl1YZzFWCPpm8LR63HjxGz/Y2Gt1xtfP1AawfbrqcmvP3sucMmlTJ0eTnaY+8Uu5RkqZ6t5e97k9kqLejzJxkLVmwbwcE+3bg/JJUJn8OjXrrnjJ/g7k/Kr9+yRDWgcqj1EJwoVAo8c/BwQFjxozBhg0bsGzZMnX3kahodeQN9F8Ti0FrY6VKnSjr0z/+kZvNlbD7KrYxv83VvDI8qW3AlZxH8AyLwva/+NW5U4a81A66TvnAJ4+XKl+qm6f1Vfg2tbTMUy4fZ/b8OJT53nCFRqQgYP0pLD30N6/2XNP+sthbGuvFYm5VqFSwtzlHR0ekpvLPh0KUl1FUgdOphazFdmOu58MzLAo/XMhCcUUNCsprVE5IJ7InQfZalO/js1S639ZgyynZgevGGO7M48p6xc8DWeGBGN/TGfZtTTG+Z2P6CH1IN8CWx6uPnN2eivhawcBJkdprrU304uGs5w1tvR+RTFyaV/oUh5Ly4BkWhUiOdWz93LlLZDXXs4NqSS/1gVJrmq5dk0x0xjAM8vPzER4ejtra2labq0kba5oeVdVg4YEUiRxMw7wdsfVVX4laTzN2JrCWpQCUn1MetfE07pZUyT3v7UhTdfJM3hbPWl6gf8d2BpUcUZ22xNzGV3/JLsb59bS+CPZVfks/1++9qWXjuimUSkOTQiNScDGrBP6dHLBxynO67g4AYPOJ23Lz5ygzgrDtVBrO3y3GUG9HzBvZBdtOpeGni9moq2/AJF+3VrFmTJdU2UGnzJqm00tG6OVubb7f30ovBBcIBGh+00GDBuGHH35A9+6GO/SmCm0ETTN3X8L59GLUN/ndGwsEGNLFAftmDZSqJcRGmS+Hw0m5rNmBAdnbtUmjHp/8iSe18ut4tDE1wq3PXtRij/RHl4+iUcdSPkGVhaNcv/em/DvZ4QBHWRdNk/c+UzV4VCdVU29cSC/C9F2XeLXdPbO/3BIvRHlc6SSm9HNjDdYPXsxWONfbnjcHYGQ3J4Vuow0aXQiemSmZIt3IyAiOjo6wsGCv9E5Uk1FUITPLdz3DiIvtKlIu4lxakcJBE1fABMjerk0acRVU7tm+de66jLiUwxowAcDSQ38rPdpiaiTAE55t9WGRt7z32cKIFL0JmlSdiuMbMAHArH1X9G63VUvAlbg0IaOY9fwrfh743+VchYrzspXqMgRKrWny8PCQ+Ofu7k4BkxZkP5Q/JQYA839SrJ6Yol8OfHPe0NScfFwFlVvr1FxCJvcaO64PcDampvw/6nQ9Ncf1PuO7UFefbTul+KaHNceUywdE5PPztGc979+JOyu7IkXiuUp1GQLeI01ff/017ztduHChUp2Rpa6uDqtWrcLPP/+MgoICtG/fHm+88QZWrFghzhfFMAw+/fRT7Ny5E6WlpfDz88P27dvRs+ezrKLV1dVYsmQJDhw4gCdPnmDUqFHYsWMH3NzcxG1KS0uxcOFCREZGAgCCg4OxdetWtGvXTm3PRxUedpas57mCquYU/XLgU05hToCnQvfZGi0d21Xmou+lY7vqoDf6gU/qgcc8CubK08/dljVXlchAT8UXt6qbqlf/huD8XcWfQ3y6YqVaCDeuxKV8R3aHdLLF+Qzu3XR8S3XpM95rmry8vPjdoUCAjIwMlTrV1Oeff44tW7bgxx9/RM+ePXHlyhW8+eabWLNmDd5//30AwPr16/H5559j79696Nq1K9asWYOzZ88iNTUV1taNu2/effddHD16FHv37oW9vT1CQ0Px8OFDJCUlwdi4sYL9iy++iLy8POzcuRMAMHfuXHh6euLo0aO8+qrrNU319fWci79FIuYMgl9n9quM5lSd/yaSXt5xHjfyyzWep8lQcK1pElF2wTGfRav6MAXUGt5n206lYVOsYrtFZwd40aJwDZCXuFTR9XN83l/6vNFFowvBtSkoKAjOzs7YvXu3+Nh//vMfWFpaYv/+/WAYBq6urli0aBGWL18OoHFUydnZGevXr8fbb7+NsrIyODo6Yv/+/Zg2bRoA4P79+3B3d0d0dDTGjRuHW7duwcfHB4mJifDz8wMAJCYmwt/fH7dv30a3bt04+6qNoKmsqhYLDiTL3T3H9sJ1sTHDTH8vlaYf1FGriBBZbt4rQ9C2ePCIm5R6rXUKiwLbUnBTIyBNTzYwtIb3maI7r1rK89ZXSw/9jYSMYqV3ao7edAbpLClwAMDUGEj7XD//jny/v9Wap0kTAgICcOrUKdy503hV8vfffyM+Ph4TJjQWc8zMzERBQQHGjh0rvo25uTmGDx+OCxcuAACSkpJQW1sr0cbV1RW9evUSt0lISIBQKBQHTEDjbkChUChu01x1dTXKy8sl/mma0NIU+2YNxOklI7DnzQE4vWQE9s0aKE43IK9m0MrAHkj8aIzK6zXk5bxRNBdOa7c68gZ8Pz0B309P0FqNf/l0ECJjXSA2vtwHJuz5LhH0FXv9S1l6ceR8ek6FzOPq1hreZxFz+Cctbcm1/fRBaEQKTt4sgADAQC/lcik95VE+zdHKXKn71idK7Z4DgLy8PERGRiInJwc1NTUS57744guVOyayfPlylJWVoXv37jA2NkZ9fT0+//xzvPrqqwCAgoICAICzs+R2VGdnZ2RnZ4vbmJmZwdbWVqqN6PYFBQVwcpLeBunk5CRu09y6devw6aefqvYElSSvEO9bQzvhraGdlK4lxEVUTkHVq5LWKuZ6Pub+dFXi2K74TOyKz6Rt1f+a0t8dHx/+B6iXP+SkTLmOdlYWAB7LPa9P0waG/j4LjUhB7K0HMDMW4K2ATjIv1vw62yMrPBA7TqfjXFoRhno7oqerDT45ch33Sp/AzESAGYM8aUpOg5qntih9Uoelv17D0l+vIXpBAHw6CGXebnXkDVzIKEZAF0fx34dPyaL1BvQalkepoOnUqVMIDg6Gl5cXUlNT0atXL2RlZYFhGDz//PNq7WBERAR++ukn/PLLL+jZsydSUlKwaNEiuLq64vXXXxe3a1owGGhcHN78WHPN28hqz3Y/H374IT744APxz+Xl5XB3149ik5qu52RIH+D6pHnA1BRtq36mi4MVbj6QHxjZWSr20SUvXUdTIbsS9a4OmqG9z2Tll9pwIhUbTqTKXUc5b2QXBPZpj8Cvz6Gi+tloxZM6BrviMzGmh7PC6y8JP2wpZIK3n0f62gkSx5pf9N0uqBBf8HEtKrcwNdKLdB6qUmp67sMPP0RoaCiuX78OCwsL/Pbbb8jNzcXw4cMxZcoUtXZw6dKlCAsLwyuvvILevXsjJCQEixcvxrp16wAALi4uACA1GlRYWCgefXJxcUFNTY1UoeHmbR48eCD1+EVFRVKjWCLm5uawsbGR+EeIPKsjuevy0VRdI65yHQWPFdtJx2dnaUreI4Xuk0hj+xKe9n2i3HOTtp+XCJj43o4ojyu1RV0Dg0NXciWOybvom7WvMd2NvOljSxMgdhH7e9pQKBU03bp1SzzKY2JigidPnqBt27ZYvXo11q9fr9YOVlVViVMLiBgbG6OhoXFJp5eXF1xcXBAbGys+X1NTg7i4OAwePBgA0K9fP5iamkq0yc/Px/Xr18Vt/P39UVZWhkuXniVcu3jxIsrKysRtCFHFBR5bxWlbdaOIS7JrGzYVsov/lylXug4A6KtHa5oMEZ88brIKhr+7/zJKOdJJqKvQOHmGTwqZpqkhuC761hy7iWDfDsgKD8SUfm5waGsGH5e22D9rIG6uCYS7Pfd70BAoNT1nZWWF6upqAI0Lqu/evSvOiVRcrN4cIhMnTsTnn3+Ojh07omfPnkhOTsYXX3yBt956C0DjlNqiRYuwdu1aeHt7w9vbG2vXroWlpSWmT58OABAKhZg1axZCQ0Nhb28POzs7LFmyBL1798bo0aMBAD169MD48eMxZ84cfPfddwAaUw4EBQXx2jlHCJfBnRxwu4B9LQ4Vi23EJ9mlIiNDnRzbYpi3I+sUnb5NzRkaPl/CsqoQ/HmDO3/W4eR7Ok862tLwWYM0pPOz5JZcF31NL/gMbVpZEUqNNA0aNEhclDcwMBChoaH4/PPP8dZbb2HQIPV+8GzduhUvv/wy5s2bhx49emDJkiV4++238dlnn4nbLFu2DIsWLcK8efPQv39/3Lt3DzExMeIcTQCwZcsWTJo0CVOnTsWQIUNgaWmJo0ePinM0AcDPP/+M3r17Y+zYsRg7diz69OmD/fv3q/X5kNZrZXBPzja06LVRThH71mVA8ZGhra/6orOj7GzE8nadEv64sksD0lUIJmyJ43XfaUos/CfsNnPsxDQxEmBK/2drdAdzZAdvLRd8SuVpysjIQEVFBfr06YOqqiosWbIE8fHx6NKlC7Zs2QIPDw9N9FXvaSNPEzFsp24+EM//N0e75/gVhBZRdtF8ZnEl3vs5CVkllXi+oy2NMKkRV+6l5n+zrh9Ho4Zll2RTLSGpp76Rl9gSgMzdcy05f1iLSW5pSChoInytOXYTvyU1LrL8Tz93GmH6F9+EhysDe+CtoZ003BuiKLYvYVm75yZsiWPdJdmUm60F4pePUrWLRIalh/7GyZsFaGthgoWjukqMMDUl76KvJVzwaTRoevPNN/Haa6/hhRde4NzW35pQ0ESI8rjKhwDA0C72NDJkAJYe+hsnbhSw5mkS4Rso00iT/lhz7Cbi04sk8jQZOr7f30otBC8pKUFgYCDs7e3xyiuvICQkBH379lW2r4QQwrmQ2M3WggImA7FxynO8A5zertb45778pKNN75Poh5YSKClDqYXgkZGRKCgowH//+18kJSWhX79+8PHxwdq1a5GVlaXmLhJCWgOuhcQ9XNjLoBDDdHThMM42Lal8DDFsalnTlJeXhwMHDuCHH35AWloa6urq1NE3g0PTc4Sohm2qxsRIIJWhmLQM8tbKDO5kh1/m+uugR0RRc3+8jKTcUgzwsMO3IYZXK1BrBXtra2tx5coVXLx4EVlZWXKzZxNCCBuu5HmyMhSTlmGUjzOywgMxO8AL3V3aYnaAF7LCAylgMgA/JWTCMywKMbcKUVJRi+M3HsAzLAoHL2brumsaoXTB3tOnT+OXX37Bb7/9hvr6ekyePBlHjx7FCy+8oM7+EUJaCT4Z0786dQcA5O7uIZoRl1qI5b9eQ0lVNZ7r0E5jxY1b81oZQ7XiD9mln8IOX8crfi0v/ZBS03Nubm4oKSnBuHHjMGPGDEycOBEWFhaa6J9Boek5QpS3OvIGfriQxautiZEAkfOHyK3CTtQju6QS47acxdO6BqlzS8d2xfwXvHXQK6Iv5v54GTG35Gd0H9/T2WCm6jQ6Pbdy5Urcv38fR44cwZQpU+QGTHl5eeIacYQQwoZPxnSRugYGwdvPa7A3BGgspCsrYAKAjTF3tNwbom+ScktZz1/OfqilnmiPUkHT3LlzYWtry9nOx8eHdtMRQnjbPZP/VSmtcdKsuNRCzkK6L++gwLU16+fOHgcM8LDTUk+0R+WF4Gwo2bjmRVzKwaKIZPryIC1C0wXBLjbmsDJl/4j6pYUuNtUHfAoi38gv13xHiN7a+foA1vOGMjWnCI0GTURz/sl7hC4fRWP57//gSPJ9LP31Grp8FI2b98p03TVCVLYiyAeJH43GyonsU3bJuWX0utcQPgWRe7antZutXfhLvRQ6bugoaDJQL+24gLoGyZE8WudBWpppAzvCxIi9VBO97jVjeDcn2FqasrbR1C46Yjhe8fNAVnggxvd0hn1bU4zv2Tha3BJ3zgEqpBwg2rPtVBr2JmShurYe43u1R38PW6mASUS0zoO2ZJOWInL+EARvPy/3NQ9o/3W/7VQafrqYjbr6BkzyddPJVvltp9Jw/m4xhno7stZ2U0Xk/ACM+TIOT2tl754jRKQlTsXJopaM4PLY2NggJSUFnTq1jmrk6k45cCG9CNN3XVL4dpN8XfHlNF+VH58QffLS9ngk58qfhtPG657tPamtSu/y+hAxZxD8OrOXolHWubQiLD/0N4oqNZuniRgWQ88C3pTWMoKzoYXgqlEmYAKAIZ0d1NwTQnTvlQEdWc9r43XP9p6UVQZEm32Y9n2ixh5zqLcjLnw0GmmfB1LARFpdFvCmNBo03bx5Ex4eLXNeU9O2nUpT6nYmRgKamiMGaXXkDYz/Mg5rjsnOMJxa8FjubbXxuufznpTXd231YcfpdI0+PiEAexbwlo73mqbJkyfzvtPff/8dAODuTl/eyjp/l7ukRHOiLMmEGJKY6/mY+9NV8c+3CyqwKz5TPN3V/Hxz2nrd83lPxqcX6bQP59KKNLa+iRCgcUqOzTv7rxj8VB0b3iNNQqFQ/M/GxganTp3ClSvPhqOTkpJw6tQpCIVU1kAd+Ew1TOnnho0v98EkX1dsfLkP0tdOoLISxODIC4hE011sAZORAFp73fN5TwZ0cdRpH4Z6a/bxCWmNWcCb4j3StGfPHvH/L1++HFOnTsW3334LY2NjAEB9fT3mzZtHNdfU5L1R3tgUy16mYOOU5wBQ8VJiuFZH3mA9P+HLs6znGxhobdccn/ekpnfRcfWBRpmIpvVzt2WtN9cSs4A3pdSaph9++AFLliwRB0wAYGxsjA8++AA//PCD2jrX2kXMGST33NfT+mqvI4RoyIUM9umm9KIKzvtQZipbWWzvSUVKwGiiD2x9I0RdWmMW8KaUytNUV1eHW7duoVu3bhLHb926RQV61civsz2ywgOx43Q6fjifKc7TJBphIsTQDe7kgNsF8gOjLo5tcZNlATig3d2iTd+T+xKydJKnqWkfzqUVaTRPkywzdibgau4jOLQ1w9rJfVSaEgyNSMGJG/kwEggwdUBHneS7IooLf6mXzEXfLTULeFNK5Wn64IMPsHfvXnz00UcYNKjx6iYxMRHh4eGYOXMmvvjiC7V31BCoO08TIeq27VQaIq7kAmDw6kAPvZjO8QyLknsuKzyQ9byoDVHc0PBTuFf2FO7t2iBu+Quc7XfGpWPtn6lSxy2MgNjQkXC3t+T92IeTcrH40DWZ57SV74qo7p39V3A5+2GrytOkVNDU0NCATZs24auvvkJ+fj4AoH379nj//fcRGhoqMW3XmlDQRPQVW1JGTSZF5OPUzQcycxyJvjzlnW9qZWAPvDW0dSTRVdVHv6bglyv3pI7PHOSO1ZP6yL0dW/Bqa2mK5JVjefeBAmGibzQaNDV/IAAUJICCJqK/DOFLas2xm4hPL0JAF0eZ0zTTv7+AC3fl79zRh+dgCLhG9mSZsTMB5zPYd0XtnzWQ11RdaEQKfkuWDtqamh3gRVN1RKu0lhHcxsaGAgRC9BifpIz6kBRxRZAPji8aLvfLki1gAoCQXZrLiN1SDA0/xXp++Pq/ZB6/dl9++RqRqznsfx+Ri1klnG00ne+KqCajqAJzfryMFzad1nhCV32jVND04MEDhISEwNXVFSYmJjA2Npb4RwjRH3x2l51L0+8vKa6EegD/L+3W7F7ZU9bzuY+eyDzex5U7D9bzHW159cHPk3sqWNP5rohyHlXVYMJX5/DC5jjE3ipERnEVdsU3llQ5dfOBrrunFUrtnnvjjTeQk5ODTz75BO3bt4dAIFB3vwghajKkswMSOKZW9D0pIldCPQDwtLfSQk8MWwehBXIfyQ+c3Nu1kXn857n+nGua+L6GNk/ryzk9R1Nz+mnhgRTczC+XeW7WviutYopcqaApPj4e586dQ9++fdXcHUKIuvFJyqgPu+jYcCXUA4BtM/ppqTeG61zYKNbgh20X3crAHlgddUvquIUREDk/QKF+fD2tLxZGpMg8p618V0QxGUUVOMsxIr3m2M0WH/AqNT3n7u4OFdePE0K0iC3xoSEkReRKqNfZ0QpeDjTSxMfMQbKzp8s7LvLW0E7ICg/E0C72aGNqBHdbC+yfNRC31wYqlG4AAIJ9OyArPBBT+rmhrZkRbMyNMTvAC1nhgZRuQE9lP6zibNMa1qIptXsuJiYGmzdvxnfffQdPT08NdMsw0e45ou92nE7HgUs50Kc8TXwdvJgtM6GeCYBXB3XEZ5N6a79TBmz4+r+Q++gJ7zxNpHXLKKrAC5vjWNsY8q5HjaYcsLW1RVVVFerq6mBpaQlTU1OJ8w8ftuyCffJQ0ESI5r2z/wrOpRWhska6+gAlRiREc2buvsQ6RWfIa5r4fn8rtabpyy+/VLZfhBCikm9D+stdl9NaFqMSogtbX/XF9F2JuHFfejF4a1mLpnJyS/IMjTQRon4ZRRXIflgFT/vGdUurI2/ghwtZctsb8hQBIYYgs7gSU749j5KKWnS0bRnTuxodaWrqyZMnqK2tlThGAQMhRFWPqmqw8ECKxHTAMG9HZJfIL/ALAH/dfkBBEyEa8lNCJlb88SyhZXbpE3iGRSH8pV54xc9Dhz3TDqV2z1VWVuK9996Dk5MT2rZtC1tbW4l/hBCiqoUHUnA+XTIx5/n0YlRU17PerouTtSa7RUir1jRgakrWJo2WSKmgadmyZfjrr7+wY8cOmJubY9euXfj000/h6uqKffv2qbuPhJBWRpQTpr7Z6oF6hkFJZQ3rbT+c0EOTXSOk1eLKzP/OfvbC2i2BUtNzR48exb59+zBixAi89dZbGDp0KLp06QIPDw/8/PPPmDFjhrr7SQhpRbhywrjbtkFuqXTJDx8Xa8rXRIiGcGXmv5zd8nfOKzXS9PDhQ3h5eQFoXL8kSjEQEBCAs2fPqq93/7p37x5ee+012Nvbw9LSEn379kVSUpL4PMMwWLVqFVxdXdGmTRuMGDECN27ckLiP6upqLFiwAA4ODrCyskJwcDDy8vIk2pSWliIkJARCoRBCoRAhISF49OiR2p8PIYSdhx17ssTc0ifwcZGchhvm7YgDc/012S1CWrV+7uzLbwZ42GmpJ7qjVNDUqVMnZGVlAQB8fHzwv//9D0DjCFS7du3U1TcAjYHMkCFDYGpqij///BM3b97E5s2bJR5nw4YN+OKLL7Bt2zZcvnwZLi4uGDNmDB4/fixus2jRIhw+fBgHDx5EfHw8KioqEBQUhPr6Z+sjpk+fjpSUFBw/fhzHjx9HSkoKQkJC1Pp8CCHcOjm2xTBvRxiz1LW8WfAYp5eMwJBO9nC2MYOzjTmElqZy2xNCVMOVmf/bkJafdkCplANbtmyBsbExFi5ciNOnTyMwMBD19fWoq6vDF198gffff19tHQwLC8P58+dx7tw5mecZhoGrqysWLVqE5cuXA2gcVXJ2dsb69evx9ttvo6ysDI6Ojti/fz+mTZsGALh//z7c3d0RHR2NcePG4datW/Dx8UFiYiL8/PwAAImJifD398ft27fRrVs3zr5SygFC1KesqhYTvj6LeywFZmX5elpfBPt20FCvCGnd5GXmN/Tdc3y/v5UaaVq8eDEWLlwIABg5ciRu3bqFAwcO4OrVq2oNmAAgMjIS/fv3x5QpU+Dk5ARfX198//334vOZmZkoKCjA2LFjxcfMzc0xfPhwXLhwAQCQlJSE2tpaiTaurq7o1auXuE1CQgKEQqE4YAKAQYMGQSgUitsQQrRHaGkKawvFl13KKwRLCFHdK34eyAoPxPiezrBva4rxPZ2RFR5o0AGTIlTO0wQAHh4e8PDQzC8sIyMD33zzDT744AN89NFHuHTpEhYuXAhzc3PMnDkTBQUFAABnZ8nSCc7OzsjOzgYAFBQUwMzMTCodgrOzs/j2BQUFcHJyknp8JycncZvmqqurUV1dLf65vFw6SyohRHmDOzngdgF7XiZZlh76GxunPKeBHumPydvicbOgHB3aWWLXGwNoATzRqtYwFSeLUiNNAHDq1CkEBQWhc+fO6NKlC4KCgnDy5El19g0A0NDQgOeffx5r166Fr68v3n77bcyZMwfffPONRDtBs7UPDMNIHWuueRtZ7dnuZ926deJF40KhEO7u7FXCCSGKWRncU6nbJWQUczcyUF/FpsIzLApX88rwtI7B3eJKjNx0BsPC/0JZVS33HRBClKZU0LRt2zaMHz8e1tbWeP/997Fw4ULY2NhgwoQJ2LZtm1o72L59e/j4SGb37dGjB3JycgAALi4uACA1GlRYWCgefXJxcUFNTQ1KS0tZ2zx48EDq8YuKiqRGsUQ+/PBDlJWVif/l5uYq8QwJIWyUqWnl38lBAz3RD1tOpcs8nvPoCRYcSNZybwhpXZQKmtatW4ctW7bgwIEDWLhwIRYuXIhffvkFW7Zswdq1a9XawSFDhiA1NVXi2J07d8TTgV5eXnBxcUFsbKz4fE1NDeLi4jB48GAAQL9+/WBqairRJj8/H9evXxe38ff3R1lZGS5duiRuc/HiRZSVlYnbNGdubg4bGxuJf4QQ9Rrl07hmYnaAF7q7tMXsAC/OorwtdWpu8rZ41vNn04qQWVyppd4Q0vootaapvLwc48ePlzo+duxY8Q42dVm8eDEGDx6MtWvXYurUqbh06RJ27tyJnTt3AmicUlu0aBHWrl0Lb29veHt7Y+3atbC0tMT06dMBAEKhELNmzUJoaCjs7e1hZ2eHJUuWoHfv3hg9ejSAxtGr8ePHY86cOfjuu+8AAHPnzkVQUBCvnXOEEM1qXk/u62l9ZS76/npaX+10SAduPXjM2SarpNIg1zc1L8xM9FtGUQXWRd/C3aIKvNDdudXUe1QqaAoODsbhw4exdOlSieN//PEHJk6cqJaOiQwYMACHDx/Ghx9+iNWrV8PLywtffvmlRNbxZcuW4cmTJ5g3bx5KS0vh5+eHmJgYWFs/S363ZcsWmJiYYOrUqXjy5AlGjRqFvXv3wtjYWNzm559/xsKFC8W77IKDg9U+3UgIUY9g3w4I9u2ApYf+RkJGMfw7ObCOMM3YmYBr98vQ160d9s8epMWeqk8PZ2tczStjbeNpb1gBh7zCzFtf9aW8W3roUVUNpn9/ETfzn218yojPxK74TOye2R+jfGQvZ2kpeOdp+vrrr8X/X15ejk2bNmHIkCHw92/MwJuYmIjz588jNDQUK1as0Exv9RzlaSJE/+yMS8faP1Oljq8M7IG3hnbSQY9U4xkWJffcMG9H7Js1UIu9Ud3M3ZdwPr1Yos6gsUCAIV0cDO65tAYzd1+SCHCb45o611d8v795B02isilcBAIBMjIy+PWyhaGgiRiKydvicevBY/Rsb4Nf5w3RdXc0ii3I0PcP+G2n0vDTxWzU1Tdgkq8bVgT5YPtfadgYc0eqbcd2bXB04VCDGp3JKKrAC5vj5J4/vWQETdXpEa6/FwDMDvAyyKk6vt/fvKfnMjMz1dIxQojufBWbKrH76krOI3iGRWHp2K6Y/4K3DnvG3+rIG7iQUYyALo6cH84zdiawng/ZlaiXU3UX0oswfdcliWO7/p0CaWsOmJsApgIj1DGMQedp4irMbKjrs1oqrr8XAMSnyx+FagnUktxSHhsbG6SkpKBTJ8MbAiekJZK3XX1jzB29Dppm7EzA1ZxHeFLXID52u6BC5jqKuNRCpOQ9wvMdbXHtPvv6n5S8R5rqskqaB0xNVfybT7cajb+LGX4dDTaw4CrMbGjrs1o6rr8XAAR0cdRCT3RH6eSWfChR1o4QoiFc29Vf3nFeSz3hb2dcOjzDonA+46FEwNTUrH1XAADZJZXwXR2D1/dcxpbYNITsvoQnNfUybyPS162durussm2n0hRqvzrqloZ6onnyCjMbCwQY5u1osMFgSyX6e7ExxKk5RWg0aCKE6A+u7eo38vWnDFBGUQVOpxbKXMAty5pjNzFp+3mUNsuILSfOEtPHqbnzdxXPZh6yK1EDPdGOra/6YkgXyWSkQ7o4YOurvjrqEWGz9VVf9HSVveZHmUS0hkaj03OEEP3BtV29Z3vdb16Qtf2cj+M38qUCJi4rA3so1F5bhnR2QELGQ4Vuo6/TjHwILU2xb9ZAZBZXIqukkvI06TmhpSmiFg5FZnElwqNvIa3wMeVpIoS0PL+/F8C6k0wfdtEtPJCi1EJSp7bmyCt9Kvf84jHeuJL5ECl5j/Q+T9N7o7yxKVZ6dxwbXU0zjtp4GtmlVfCyt0Js6AiV7svLgYIlQ+LlYIXvWsHIUnManZ7jKphLCNGupWO7KnRcmzKKKnA2rQgNSiyFfH80e/+f72iL/bMH4Z9V4/U6YBKJmKNYH7X9nD794x94hkXhbkkV6hqAtKJKeIZF4fNjN7TaD0K0jRaCE9KKzH/BG1nhgejfsR3amBqhf8d2yAoP1Iudc3y2M8uye2Z/DO/mBGM512i2lqYYyrF4Vd/4dbZHVngglo3rBhcbc7RrYyy3rS6mGfck5Mg8/n18lnY7QoiW8U5uqYz4+HgMGDAA5ubmmnoIvULJLQlRHp/Eec3tntkf+WVVWPHHTZnnbS1NETk/AO723FulDcGhK7nYcPw2yp/WYKCnvU5GzUZtPI27JfIDXG9H1afqRFpC6RtiGNSe3PKDDz7g/eBffPEFACAgIID3bQghrVsnx7bo0M4C9x7JX5vUnCjdgDylVbVIziltEUGTZBb3MTrrR3Yp+4hgZkmlyo/RvPTNufQSeIZFGWzpG9Jy8A6akpOTJX5OSkpCfX09unXrBgC4c+cOjI2N0a9fP/X2kBDSakQvHIbnVseo9T4XRqTAt6OtwQZO+pbF3cPWknWkyUsNCSnlpZpYHXWLgiaiU0pNz33xxRc4c+YMfvzxR9ja2gIASktL8eabb2Lo0KEIDQ1Ve0cNAU3PEaK6UzcfcI4gKcpYALzQ3Ql3iyoMbnu0PtbO01SfQiNScPTafdTUy/9aGtpFN9OSpGXj+/2t1ELwzZs3Y926deKACQBsbW2xZs0abN68WZm7JIQQAMAoH2d8Pa2vWu+zngFibxUio7gKu+Iz4RkWhVM3H6j1MTRBX7O4zwnwVOg4l8NJufAMi8JvyfdYAybAsHNSEcOnVNBUXl6OBw+kP3AKCwvx+DF71mFCCOES7NsBWeGBmNLPDW1MNbPJV92jWZqgr1ncPw7qiazwQHg7WsHEqHHxd1Z4ID4O6qnU/S0+dI13W30sfUNaD6U+jV566SW8+eab+PXXX5GXl4e8vDz8+uuvmDVrFiZPnqzuPhJCWqmNU57DqonKfRHzseaY7F13+qKHszXreV1ncY8NHYH0tYEq7ZYLjUhRqD1NzRFdUipo+vbbbxEYGIjXXnsNHh4e8PDwwIwZM/Diiy9ix44d6u4jIaQVmzawI0yMNJMoV5ns49r0+3vsO5D1IYu7qi5mlfBuq6+lb0jroVTQZGlpiR07dqCkpATJycm4evUqHj58iB07dsDKitLgE0LUK3L+EGgibAroov9JL3WZxX3GzgT0WvknJnx5FpnFqqcSkMXP0571vLGgcfF3Vngg7ZwjOqdScsv09HTcvXsXw4YNQ5s2bcAwTKsunUK75wjRnLjUQry+57Ja71NXu8+U8fKO87iRX/5vnibNjjA1z5Mk4mVniSPvBUBoaarWx9PHHYKkddHo7rmSkhKMGjUKXbt2xYQJE5Cfnw8AmD17dqtNN0AI0azh3Zxgq8Yv690GVmz013lDcOuzF7UyJScvT1LmwyosOJAs85wq5O2WVPcuSkJUpVTQtHjxYpiamiInJweWls8Sxk2bNg3Hjx9XW+cIIaSpyPkBSgVOY3o4YZyPMzo5WGJ2gBeywgMxysdZAz00fDN2JrCeP5tWpPapuqa7Jd1sLTClnxuywgMR7NtBrY9DiKp4ZwRvKiYmBidOnICbm5vEcW9vb2RnZ6ulY4QQ0py7vSWSV47FubQihOy+xPt2Y3u6YEp/dw32THO2nUrD+bvFGOrtiHkju2j88a7dL+Nsk1VSCS8H9a9f3TjlObXfJ1G/iEs5SMgswZDODgb7vlKWUkFTZWWlxAiTSHFxcaspzksI0Z2h3o6ImDMI075P5NXeED/YL6QXYfquZ4FhQsZDbDiRiog5g+DXmX3xtCr6uApxPuMhaxtPNZRKIYbnn7xHeGnHBdQ1NC6FPpJ8Hx/+/g8i5w+BTwehjnunHUpNzw0bNgz79u0T/ywQCNDQ0ICNGzdi5MiRauscIYTI49e5cUdVUG8X1nbRCwyzcHjTgKkpvoGisn6e6896fpi3o0ZGmYj+axowidQ1MAjerpvM9LqgVNC0adMmfPfdd3jxxRdRU1ODZcuWoVevXjh79izWr1+v7j4SQohc22b0k7uoe/fM/gZ5BbztVBrr+R2n01nPq0pePiQvO0tsfdVXo49N9FPEpRypgEmkroHBoSu5Wu6RbigcNNXW1mLevHmIjIzEwIEDMWbMGFRWVmLy5MlITk5G586dNdFPQgiRa5SPM7LCAzE7wAvdXdoa/GLv83eLWc+fS9NsUs63hnZCVngghnaxR1szI/i4WOP0khE4vWyk2tMNEMOQkMmehJTrNdtSKLymydTUFNevX4e9vT0+/fRTTfSJEEKUsiLIR9ddUIshnR2QwLKuaKi3dpJyUskSIuLvZY8jyfflnh/S2UGLvdEdpabnZs6cid27d6u7L4QQQgC8N8qb9bw2dtER0hRbOSMTI4FBbrZQhlK752pqarBr1y7Exsaif//+UqVTvvjiC7V0jhCiOT1WRONJHQNLEwFurpmg6+6QZuTtDoyYY3ijP6sjb+BCRjECuji2mNHA1ihy/hAEbz8vsbbJxEiAyPmGXwORL6XKqLDtkBMIBPjrr79U6pShao1lVGbsTMC1+2Xo69aOhvINxLRv4nExWzoXT0BnW/w0Z7AOekTY7DidjnNpRVrL0yRPxKUcfHf2Lqpq6hDUpwOv4Cfmej7m/nRV6vjumf0Ndr0ZAQ5dycX5u8UtKk8T3+9vlWrPEUmtKWiSV5tqZWAPKqqp51pKna/QiBRczCqBfycHSoqoQf/kPcL/bT8PWRunuIKflvJaIy2fRmvPESKvNtXqqFta7glRRI8V0aznfTjOyxKXWojFEcn4ICJZ47u6AOBwUi48w6LwW/I95JU+xaGkPHiGRSEy+Z7c26yOvIHxX8ZhzbGbGu9fS/PSjgsyAyYAmLXvitzbrY68wXq/9LcghkipNU1E+0ZtPI27JVUwAvCffm46vbLmqk0VsiuRpur01JM69oHlKo7zTWWXVGLi1niUP60TH/s9+T5sLIwRtWAY3O2lqwYoYtupNBz75z7srczwzogu4h1jiw9dk9l+YUSKVK2y5tNDtwsqsCs+k6aHeGLLzSOy5thNmVN1FzLYt6DHp2s+wCZE3WikSc99+sc/8AyLwt2SKgBAA8DrylqTuGpTpeQ90k5HiMLamMje/SJiyXG+qUnbz0sETCLlT+sRvD1e4b6JXEgvgmdYFDbF3sHtggqcv/sQIbsvodfKaLzNMrIBAEsP/S3xs6z1NAD7CAl5his3DyA/+BnciX0LekAX7aRNIESdKGjSc3sScuSeWxiRor2ONNHHlT3Dcl+3dtrpCFHYLY5dcnx30cWlFqK0qlbu+dKqWqWn6uSVD6moYRBz8wHrbROajG7Q9JDq/L24a9zJC35WBvdkvR3toiOGiIImPTZq42nONs2vrLWBqzYVTc3pt4DOtgodl4XPaOLVnFLe9yfCVT6Ea/LQv8noBk0PqY4tN48IW/DDVt6GGK7QiBQErD+lk+8fXaOgSY9ll1Zxtkng+GLQFHm1qeQdJ/rjpzmDkRUeKJ6KszQRICs8UKF0A3xGE5/vyD8IE1G1FEPTtX40PaQekfOHQF7cxBX8tLTyNq2dMpswWhpKOaBG6k45IFr8zWaKCovCJ2yJQ3pxJbo6tcWx94cpdR8huxKRkveI8jS1Qr6rY+RO0dlamiJ55ViF73PbqTRsir3D2mbe8M7YEXdX6vjX0/pKLQSnLe/qc+hKLr49cxeVNbW88zSRlqUlv58oT5MOaCJPE9uLFFDuhbo++ia+OZspdXzByM4IHddd4fsjrVNuSRUCt56TWgyu6u45tte8tbkJ/vl0HIDGqemEjGLWPE2nbj6Queibds/pv7k/XsblrIfwcmiLzdP6wsvBivtGRGNCI1LwG8uIkioX8PqgReZpWrduHQQCARYtWiQ+xjAMVq1aBVdXV7Rp0wYjRozAjRuSC0Crq6uxYMECODg4wMrKCsHBwcjLy5NoU1paipCQEAiFQgiFQoSEhODRo0daeFbs5gR4yj339bS+St2nrIAJALaelr56J0Qed3tLXFs1DvtnDcRk3w6Y7OuK/bMG4tqq8SqlG2ArE/K4+lmAtnHKc4hfPor1g5qmhwzPTwmZ8AyLQsytQpQ+qcPV3EcYuekMRm86gzKWzQdEsy5mse+k1NVSEW0zmKDp8uXL2LlzJ/r06SNxfMOGDfjiiy+wbds2XL58GS4uLhgzZgweP34sbrNo0SIcPnwYBw8eRHx8PCoqKhAUFIT6+npxm+nTpyMlJQXHjx/H8ePHkZKSgpCQEK09P3k+DuqJrPBAeDs2XmUZoTGizwoPlJqK4GPCljjW80FfnVWmm6QVG+rtiC+m9cUX03zFuZRUcTnrIev5HafTxf+fUVSBTyNvYPaPl3HoSq7c26wI8sHxRcNpSskArPhD9q7G9OJKLDiQrOXeEBFvR2vW8/4cawhbCoNIbllRUYEZM2bg+++/x5o1a8THGYbBl19+iY8//hiTJ08GAPz4449wdnbGL7/8grfffhtlZWXYvXs39u/fj9GjRwMAfvrpJ7i7u+PkyZMYN24cbt26hePHjyMxMRF+fn4AgO+//x7+/v5ITU1Ft27dtP+km4kNHaGW+0kvrmQ9f6ewQi2PQ4iyuBaDn0srwnS/jnhz72Uk5zwSHz95qxBhv13D0fcC4NOBPS0G0U9zf7zMev5sWhEyiytpqk4HzqWzvy8NeWpOEQYx0jR//nwEBgaKgx6RzMxMFBQUYOzYZwtOzc3NMXz4cFy4cAEAkJSUhNraWok2rq6u6NWrl7hNQkIChEKhOGACgEGDBkEoFIrbyFJdXY3y8nKJf/quC8eHTVentlrqCSGyDenMfsU61NsRCw+kSARMIvUMELz9vIZ6RjQtKZc7TUVWCfuFH1E/PpnhMzkuyFsKvQ+aDh48iKtXr2LdunVS5woKCgAAzs6S6xOcnZ3F5woKCmBmZgZbW1vWNk5OTlL37+TkJG4jy7p168RroIRCIdzd9b/ac/Ti4aznld1FR4i6vDfKm/X8+F4uOMuSOLOugWGdqiP6q587d5oKT3saZdI2PpnhW0swq9dBU25uLt5//3389NNPsLCwkNtOIJBMIsIwjNSx5pq3kdWe634+/PBDlJWVif/l5hrGB/WCkZ0VOk6ItslbDB4xZxCyH3LnL1M13xORbfK2ePT45E+8vEMzo3k7Xx/Aen6YtyNNzekAn8zwrSWY1eugKSkpCYWFhejXrx9MTExgYmKCuLg4fP311zAxMRGPMDUfDSosLBSfc3FxQU1NDUpLS1nbPHggXZ6hqKhIahSrKXNzc9jY2Ej8MwSh47ojKzwQvdpbw8xYgF7trZEVHkjpBoje8Otsj6zwQCwb1w3+neywbFw3ZIUHwq+zPTzsuHfmcU3xEcV8FZsKz7AoXM0rw5PaBlzJeQTPsChs/4s9g7sywl/qJfN4FwcrbH3VV+2PR7hxZYZvTcGsXudpevz4MbKzsyWOvfnmm+jevTuWL1+Onj17wtXVFYsXL8ayZcsAADU1NXBycsL69evFC8EdHR3x008/YerUqQCA/Px8uLm5ITo6WrwQ3MfHBxcvXsTAgQMBABcvXsSgQYNw+/Zt3gvBNZGnqalBn8fiweMatLcxx4WPRnPfgBAZVkfeQPT1fFiZGeOdEV0wpb/+Tys3N3P3JblTdCZGAqSv5VdDj3D7KjYVW06lyz2vqaSG7+y/gosZJZSnSU/cvFeG4O3npdY2+boJsfctPwgtTXXUM/VoscktR4wYgb59++LLL78EAKxfvx7r1q3Dnj174O3tjbVr1+LMmTNITU2FtXXjFsl3330Xx44dw969e2FnZ4clS5agpKQESUlJMDY2BgC8+OKLuH//Pr777jsAwNy5c+Hh4YGjR4/y7pumgqbQg1fxW0q+1PFp/Ttg/ct91fY4pGWLuZ6PuT9dlXkueoFh7Tgrq6rFm3sv4WqzxeDGAtDuOTXjSrDbv2M7/DpviJZ6Q3Tt0JVcxNwogLudJUL8PVtMMMv3+9sgUg6wWbZsGZ48eYJ58+ahtLQUfn5+iImJEQdMALBlyxaYmJhg6tSpePLkCUaNGoW9e/eKAyYA+Pnnn7Fw4ULxLrvg4GBs27ZN689HFlkBEwBEXLlHQRPhTV7ABDTuODOk0RmhpSl+nzcEmcWV+CkhC9kPqzCup4tBjprps8nb4jnb3MjX/13DRH2m9Hdv1e8zgxtp0meaGGka9HksCh7XyD3vSlN1hIfVkTfww4Us1jYbX+7Tqj8MibQen/yJJ7UNrG1opIm0BC2yjEpr9IAlYAKA/PJqDPo8Fl5hURi89qSWekUMzQUeJQ5oxxlproczexZoABQwkVaFgiY952xtxnqeAVDwuAYMgPvl1fAMi8LyX1O00TViQAbzKHFAO85Ic7+/F8B6funYrlrqCSH6gYImPZf48RiFbxNxRX4latI6rQzuyXrexEhAU3NEJnmB0dKxXTH/BfZEpIS0NBQ0GYBp/RUvzEtTdaS53TP7yz0XOZ+mWIhs81/wRlZ4IPp3bIc2pkbo37EdssIDKWAirRItBFcjTeZp4tr225wAQKaG8qcQw7bm2E0cu3bfoPM0EUKIOrWalAOtwYQv4xS+TXsbcw30hLQEK4J8sCLIR9fdIIQQg0NBkwG4WVCh8G0oDQEhRJ7J2+Jxs6AcHdpZYtcbA5ROUDj3x8s4nVoIBsDoHs74NkT+FDAhLQGtadJzcamFCt9GmTVQhJCWr2kNuad1DO4WV2LkpjMYFv4Xyqpqed/PTwmZ8AyLQsytQtQ2AHUNwPEbD+AZFoWDF7O574AQA0VBk55LyXvE2cbVxhyCf/+bFR5IWcIJITLJqyGX8+gJFhxI5n0/K/64Kfdc2OHrCveLEENB03N6rq9bO9bzU57vgI1T+2qlL6TlCI1IwZm0Qji3NceHgT4Y6u2o6y4pLOJSDvZeyEJNXT3+088d80Z20XWX9BpXSZSzaUXILK7knKqb++Nlzsd6Z/8Vzqm6uNRCHEm5BwGAl553M8jXYEsVGpGC83eL0dmxLda81LvF1JdTBwqa9Nzwbk6wtTRFqYyhc1tLUwqYiEIOJ+Vi8aFr4p9LKmoRsvsSLE2AE4tHwt3eUoe94+efvEeYtP086pvs+91wIhUbTqQiYs4g+HW2113n9NitB48522SVcAdNSbmlnPdzOfuh3HPZJZWYuDUe5U/rxMd+T74PGwtjRC0YZhCvwZaq+edDQXk1Rm46g25ObfG/dwZDaGmqw97pB5qeMwCR8wNg2+zFaixoPE6IIpp+IDZVVQcEb+cuzqoPXtpxQSJgamra94na7YwB4VMSxdOee0Shn7stZ5sBHnZyz03afl4iYBIpf1pvMK/Blkre50NqYYVC07ctGQVNBsDd3hKv9HeTOFbPAEM3nsbmE7d11CtiaEIjUljPl1bV4lxakXY6I8PcHy+j35oYvLP/itw2EZdyUNfAnlpux2nZ63ZaO66SKMO8HXlNw+x8fQBnG3lTc3GphTJHzUV0/RpsrbJLKtHlQ/ZcgKLp29aOgiYD8c3ZTJnHt56+q+WeEEN1MauEs83VHO6pF3VruhOrpKKWdRdWQib3c6AvXfnklUTp2K4Ntr7qy/t+wl/qpdQ5PhtbdPEabM2ySyoxfOMZ1PFIc51VQkETBU0GYMIW9uSWQV+d1VJPiCHz8+Re6/N8R+6pF3WTtxNL1i4sfy/u59DSFhRnFFXgdGqhWq7ym5ZEsTARoLODFU4vGYGzYS8otF7lFT8PZIUHYnxPZ5gaASZGwPiezsgKD8Qrfh5yb8e1sQXQzWuwNQv8+hzvtnymb1s6WghuANI5PizvFCqe/JK0Ppun9cVvyfKLOdtammo94ODaidV8F9a0gR3x8ZHrrFN0LWUX3aOqGiw8kIKzTUbOhnk7YuurviovyP11nnpqDSqazJJtYwugm9dgaxaXWoiK6npebflO37Z0NNJkALpwvFC7OrUV//+MnQnoviIa/T+LwaEruZruGjEgX8Wmyj1naaKbjQVcO7Fk7cKKnD8ExgLZ7SPmDFJHt/TCwgMpOJ9eLHHsfHqxwS/IjZwfABsL6et1Gwtj2tyiZXymSwGgm1NbhaZvWzIq2KtGuirYmxUeiJ1x6Vj7p+wvxegFAfDpIFRrf4jh4XoN6cLcHy8j5pb8rPfje8ovzXHoSi72xGeiugXmacooqsALm+VPy59eMkJvrvq3nUrD+bvFGOrtqNDf4FxaEQ5fvQeAoTxNOhKXWojX97CP9u6fNbBV/G34fn9T0KRGmgyaNp+4LXPR94KRnRE6rjvrF6KJkQDpayeotT/EsEzeFo+reWVyz/fv2E5tUzaK0sdgTtdOpxbiTZYvsz1vDsDIbk5a7JG0C+lFmL7rktRxypVlWHxXx8idLh3oaYf/veOv5R7pBt/vb5qeMxCh47ojKzwQvdpbw8xYgF7trZEVHojQcd0xY2cC623rGhiaqmvluBIb3sgv11JPpMnbbcW2C6ul87BjT/CoDwtyZQVMAOXKMjSR8wPQro30GrkBHrb4fqb0KO+2U2l4dWdCq03tQUGTgTn2/jDc+XwCjr0/THzs2n35Iwgi5+8Wc7YhLRdXYsOe7dU7MqqIpjux7Nua8tqF1dJ1cmyLYd6OMBZILt4yFgj0YkHutlNprOdb6xeqIXK3t8SG//SWOn45uxRXsp6tKbyQXgTPsChsir2DhIyH2HCisfjzxbvcaUBaEgqaWoA+rtzrlYZ0dtBCT4i+4kpsqKupuaa+DemPpBVjFd6R1VJtfdUXQ7pIvm+HdHHQiwW5XBdhG06k4tTNB1rqDVHV3J+uyjw+a9+zRLM0stiIgqYW4Oe57HPOJkYCTOnvrqXeEH0lL7GhvONEt4SWptg3ayBOLxmBPW8OwOklI7Bv1kC9qP/F5yKs6Rcu0V+rI2+wnl9z7CaNLDZBQVMLsTKwh9xzdQ0MfjiXocXeEH3UNLFhG1Mj9O/YDlnhgZj/greuu0ZYeDlYYWQ3J51PyTX13ih+r5k1x2QnLiX640IG+6hhfHoR58hia8rCT0FTC/HW0E6sO41WR93SYm+IPvt13hDc+uxFvZiSI4aLT06s+PTW82VqqAZ3Yh81rK1rgI05++hma0hJIEJBUwvCtYsuZFfrmnsmpKWY++Nl+H56ApO3n9eboql+ne3x1mBP1jYBXVrPl6mhWhnck/X83eIqnLjFvj6tJeVI40JBkwGasCUOXT+Olqo5x7WLjm/2V0KIfmhazLj0SR2u5j7CyE1nMHrTGZTJya2jTVxfuCuCfLTUE6KK3TJSC/DVkrLw80FBkwFZH30TnmFRuPmgAjX1DK7nP4ZnWBQ2n7gNgHsXHZ9imYQQ/SGvmHF6caXelFOR94Wryhcx0a5RPo1pPmYHeHG2Hd/TGf6d7LBsXDdkhQe2ukSmFDQZkG/OZso8LsoUzrWLbv/s1nVFQIghe+179un2s2lFejFV1/QLt7tLW8wO8EJWeCBG+TjrumtEQd5N6pjKY2FmjANz/VvVlFxTFDQZiAlb5NehAiCeqpO3i45tdx0hRP+cvytdrLi5rBLdB00iK4J8cHzRcJqSM2Dr/mRPPwBQzj/pUtNEL6VzXFHeKawA0LiL7q2hnRCyKxEpeY/Q160djTARYmDiUgvBpyioPpRTIS3HoycNnG1ae84/CpoMRBcHK9x8UCH3fNdmw6oUKBFiuPhs2ujlaqNXuZuIYeuxIpqzzXNuuiu3pC9oes5ARC8eznq+aS06QlqSHiui4RkWBR8eH+otBZ9NGy42ZprvCGk1ntRxj23W8GjT0lHQZEAWjOys0HFCtEFTVc+nfRMPz7Ao8Yd5VR0Dz7AovPb9BQCN5R/GfxnXIrNOD+/mBAFHm9ssI8+EKKqNCdcrDpj4nKsWeqLfBAzDUOioJuXl5RAKhSgrK4ONjeaGMYO+Oos7hRXo6tSWRpiIzlxIL5JZxDNiziC1bEP2DItSqP3umf1b1I6tt/ddwQmWordT+rlh45TntNgj0tJxvefYqk4YOr7f3xQ0qZG2giZC9AHbB6yqH649VkTzmi5Q9+PqG03+jglp7rXvLyD+bqnc87aWpoicHwB3e0st9ko7+H5/0/QcIURhmq56rkzABLS8ArFfT+ur0HFCVPHTnMHICg+EpZyputKqWgRvj9dyr/QLBU2EEIVpuuo5n/UVsrS0ArHBvh2QFR6IKf3c4GZrgSn93JAVHohg3w667hppwb4JkZ/NvbSqVuX3tyHT+6Bp3bp1GDBgAKytreHk5IRJkyYhNTVVog3DMFi1ahVcXV3Rpk0bjBgxAjduSCbpqq6uxoIFC+Dg4AArKysEBwcjLy9Pok1paSlCQkIgFAohFAoREhKCR48eafop8jb3x8votfJPDFobi0NXcnXdHdKKcSW4U7Xq+a01E5S6XUstELtxynOIXz6K1jARreBKeXE1R/4UXkun90FTXFwc5s+fj8TERMTGxqKurg5jx45FZeWzZI8bNmzAF198gW3btuHy5ctwcXHBmDFj8PjxY3GbRYsW4fDhwzh48CDi4+NRUVGBoKAg1NfXi9tMnz4dKSkpOH78OI4fP46UlBSEhIRo9fnK0rRoZ0VNAwrKa7D012uNdejusRfpJUQT3hvlzXpeHSUWAjrbKnwbykZNiOq4Ul4831Hx92ZLYXALwYuKiuDk5IS4uDgMGzYMDMPA1dUVixYtwvLlywE0jio5Oztj/fr1ePvtt1FWVgZHR0fs378f06ZNAwDcv38f7u7uiI6Oxrhx43Dr1i34+PggMTERfn5+AIDExET4+/vj9u3b6NatG2ffNLUQnG0xqImRAOlrlbsqJ0QVF++WYNr3iVLH1bV7TsRnRTSq6hhYmghwc80EnLr5ALP2XZFq19J2zxGiS76rY1BaVSt13NbSFMkrx+qgR5rVYheCl5U1jqzY2dkBADIzM1FQUICxY5/9Ec3NzTF8+HBcuNCYzyUpKQm1tbUSbVxdXdGrVy9xm4SEBAiFQnHABACDBg2CUCgUt9GFuT9eZj1f18DQVB3RCb/O9sgKD8Sycd00WvX85poJyAoPxM1/p+yoQCwhmhc5PwC2lqYSx0S751ozgyqjwjAMPvjgAwQEBKBXr14AgIKCAgCAs7PkB6azszOys7PFbczMzGBrayvVRnT7goICODk5ST2mk5OTuE1z1dXVqK6uFv9cXl6u5DOTLymXe+74/N3iVl8PiOjOvJFddFLxnKbiCNEcd3tLJK8ci3NpRbiaU4rnO9qqvFaxJTCokab33nsP165dw4EDB6TOCQSSu20YhpE61lzzNrLas93PunXrxIvGhUIh3N3VH7j0c+eeO07KeoiXtsfTiBPRuoyiCnwaeQOzf7xMrz+isNCIFASsP4Wlh/7WdVeIHEO9HfH+qK4UMP3LYIKmBQsWIDIyEqdPn4abm5v4uIuLCwBIjQYVFhaKR59cXFxQU1OD0tJS1jYPHkhn3y0qKpIaxRL58MMPUVZWJv6Xm6v+L42drw/gbJNb+hTJuWVY+us1eGlpcXjEpRzM/vEyVh+9gcziSu4bkBblUVUNXtpxHi9sjsOeC1k4easQS3+9hs4f0uYEwu1wUi48w6LwW/I95JU+xaGkPHiGRSEy+Z6uu0YIK70PmhiGwXvvvYfff/8df/31F7y8vCTOe3l5wcXFBbGxseJjNTU1iIuLw+DBgwEA/fr1g6mpqUSb/Px8XL9+XdzG398fZWVluHTpWVmIixcvoqysTNymOXNzc9jY2Ej804Twl3rxbssACN5+XiP9AIB/8h6hy0fRWP77Pzh5qxA/nM/CyE1n8NK2eJTJWDRIWqaFB1KQnPNI6ng9o9nXH9E/EZdy8OrOBMz98TLv/D2LD12TeXxhRIoae0aI+un97rl58+bhl19+wR9//CGxg00oFKJNmzYAgPXr12PdunXYs2cPvL29sXbtWpw5cwapqamwtrYGALz77rs4duwY9u7dCzs7OyxZsgQlJSVISkqCsbExAODFF1/E/fv38d133wEA5s6dCw8PDxw9epRXXzVdRuWd/VcQn1YEI4EA5dX1rG03vtxHI+ucunwUjboG2S+ZYd6O2DdroNofk+iXjKIKvLA5jrWNpl5/RH/8k/cIk3acR32D5HFrcyNELxwut9RGaEQKfmMZUaKaevojLrUQKXmPWsV6Jr7f33q/EPybb74BAIwYMULi+J49e/DGG28AAJYtW4YnT55g3rx5KC0thZ+fH2JiYsQBEwBs2bIFJiYmmDp1Kp48eYJRo0Zh79694oAJAH7++WcsXLhQvMsuODgY27Zt0+wTVMC3/2ZpXRSRjCPJ91nbamJxeMSlHLkBEwCcTStCZnElvBys1Pq4RL9kP6zibEObE1q+l3ZckAqYAOBxdQOCt8fL3ZZ+MauE9X4TMtizzRPNyy6pxKTt5yVSDrTkunOK0Pugic9AmEAgwKpVq7Bq1Sq5bSwsLLB161Zs3bpVbhs7Ozv89NNPynRTq/y97DmDJq6MzcpIyGT/sAOArBIKmlo6DzvuD01NvP6I/uC6gBKV2pA1OuHnaY+8UvkjTf6d6LWja80DJuBZ3bmWmKNJEXq/polImzawI0yM2HcGrj52AyG7pBMPqsLfizv/jqc9BUwtXSfHthjGMlRvYiSgUaYWjs8FlLxSG5s5ig3T1JxuxaUWykxqCVDdOYCCJoMVOX8IjFn+eo+f1uNcegk8w6Lww7kMtTwmV7D2nJuQRplaia2v+uL5ju2kjhsLGl+bLd3cHy+j35oYvLNfOjN5a1BdU8fZhq3UxtdyAid5x4n2UN05dnq/ENyQaHohuMjkbfG49eAxera3wbSBHfHLxcYknsm58rd6Z4UHquWxb94rw4St8Rp/HKJ53T6OQnV945XTj7MGKrXQM7O4Ej8lZCH7YRXG9XRp8SNMPyVkYsUfN6WOh7/UC6/4eeigR7rBVtqpKa6SOksP/Y2EjGL4d3KgESY90GNFNJ7UsYcE+5X8rNB3fL+/KWhSI00HTV/FpmLLqXSp40vHdsWF9GKcz3go97ZDu9hj/+xBKvdh26k0bIq9I/f8snHddJIdmvD3n21nkZT3WOa5c0tHtvqFnmzYgoXWcsEwY2cC62dNc63l92LIpn0Tj4vZ3PnVWmrdOaAF155rzWQFTACwMeYOrt1nf8FzDbnyFcGR9fnApRy1PA7RHHkBEwAEb5c/itjacdWBbC1TdVyfNc3tOC37c4voD74BU2uvOwdQ0GQwJm9j/zJrqGcfMOzr1k4t/eAemKSBS33W7WP2aRVa6CkfVx3Iy9n8R18MWR9XoULt6fWk33qsiGY9b4zGKbnklWNpFBoUNBmMWw/kjw4AgIx0KRLUMTUHAK8M6Mh6/tWBrWddhyHiyIkKgBZ6ysNVB3KAh52WeqJbP8/1V6h9S1z/0pJwrWGqB/0Nm6KgyUD0cLZmPd+zvQ1WBvaQeU7ecWW8N8qb9TytZ9Jv5sbcbdh2PbVmXHUgRclnWwNFPlPoM0G/tTFhT19jyXG+taGgyUD8/h77XPKv84bgraGdkBUeCG8HKxgLABcbc5xeMgJvDe2k1r5EzJE9aiXvONEfqZ9zL8qNTGFPnNqayasDqUh9yJZA9FkztIs9rC2M0aeD7IWz9Jmg/26tmcB6/ibH+dZG7zOCk2eWju2KjTHSO9eWju0KoLEW1Es7Logz9RaUV2PkpjPwdRNi71t+EFqaqqUffp3tkRUeiB2n08VZf+lq0nAYgX06l8pYyPeKnwde8fPAO/uv4HL2QwzwsGtVI0zNNZ/2p88EwxTQ2Rbxd6Wn5QM606hzc5RyQI20lafp5R3ncSO/HD3b2+DXec8SCVIxXcIHFUwlhMjisyIaVXUMLE0ErW6EifI06YC2giZZIi7lYPnv/7C2Ob1kBGXsJgAo3xAhhDRFeZpaGb7FdAkBVCtjMffHy+jxyZ/w/fQE5eAhhLQqFDS1EFRMlygi2LcDssIDMaWfG9xsLTClnxuywgMR7NtB7m1+SsiEZ1gUYm4V4kltA0qf1GHDiVR4hkXh4l3uoJ0QQgwdTc+pkS6n5wBa00Q0i6veGE3rEUIMFU3PtUKR84fAxEg6p4avmxBbX/XVQY9IS8FVQgSgchmEkJaPUg60ID4dhEhfOwGHruQi5kYB3O0sEeLvSYu/icq4SogAjeUyaJs5IaQlo6CpBZrS3x1T+rvruhukBennbouYW4WsbajUAiGkpaPpOUIIJ64SIgCVyyCEtHwUNBFCeLHkGJce+FmMdjpCCCE6QtNzBsL30xMofVIHuzYmuPrfcbruDmmFqurYzxdW1mqnI4QQoiM00qTn3t57CZ5hUSh90viN9fBJHTzDojD/pys67hlpbRysuGsX/n979x5VZZ3uAfy73cBmK+IEqSFsL1ASDIZml0HQ1E5YaWLLEBtEB7XWTIzjNHk6uo4zXo45lomMFU6pCw/GqEOoywujeNZohCheAi9goJgXxLwrogUCz/mjJHcIvNveywa+n7X2H7z79+79fR9e4OH3XvaibV/pkISIyBhsmpzctq8u3nP5liPndU5Cbd3+P0c2O+aDHaU6JCEiMgabJifWb862Jp9/vJnnidT2cl+fZseM+Fu2DkmIiPTHpsmJ3Tkk15grzTxPpLbFYx+Hm7nhDVTvVnKhUqc0RET6YtPkxB6wNn2evlczzxNp4eFmbpbau4uHTkmIiPTFpsmJ5TdzlRyvoiMjZL75TJPPb546SKckRET6YtPk5IaHdHVoOZEepgwJcGi51oJmZqLn9C0InplpyPsTUdtgEhExOkRrofRTku/H43O24Qrv00ROZsTfslFyoRK9u3gYMsMUszQHeaeuN1geEfAAPn1tgO55iKhlUvr3myfFtBBslMgZGX0o7l4NEwDklDb/AcOkrbkbC5F74hIiHu6MmSOCjY5DpArONKlIy5kmR7y1tgB5Jy8jzP9BLIwONSwHkZaCZmbi25rGf321dzGhaN6LOiYiAMg6cg6vf/plg+Urxj+BZ4N5WgE5J6V/v3lOUyuy/sAZ9Jy+BRn5Z1F29TukHyhDz+lbsDH/rNHRiFTXVMMEALeaeZ7Ud+Ji5T0bJgCYlMpPMWgJ1u49jZC//Au9pm/Bc4t2Gh3H6bBpakXeTD90z+V/WFugbxAiHVhdmr5fVPtmnif1XLtVjfEr9mLoos+bHDdvc5FOichRh8uuodf0LfivdYdRWV0HAXDs4k30nL4F72wuNDqe02DT1Eq81Uxj9J/pB/UJQqSTo80ceuOhOf38YXUBdh2/1Oy4nOP3/lgoMt7LyblobG52Wc5JPaM4NTZNrUTeyctNPr/7RPO/0IhamoiABxxaTuo7cbES2ccuolbB6bERD3fWIRE5au3e06ipa/r7x0N13+PVc63E0z29UXa18XOXwvwf1DENkT7u3FYgeGYmbtUIT/42wKkrtxSP5VV0zmn3103/0w0AX1++qUMS58emqZVYFNMXGU2c8M2r6Kg1Y6NknB5e7RWNWzH+CY2T0P0K6+WNDfnlTY7p5d30xye1FTw814osienr0HIiop/Lv7MHBj3SGWaT/Yn3JgAd3MyYHNELJxcM5+0GnFjMU93h0q7pCye2vzVYnzBOjk1TKzKyny9OLhiO6P5+8HvAHdH9/XBywXCM7OdrdDQiasU+eLUfwh+2PwVg4COdkTv9WR6SayE2JoSjsbbptYieekZxary55U8kJydj4cKFOHfuHH75y18iKSkJAwcOVLSus9zckkhrv3pnO765UV3/tYdrOxz5nxcMTER6i/1kNw6VX0dfv19gTlQIXlu5D6WXmz6/6eSC4Tqlo/vxefEFTEjZV//1I507tJkZJqV/v9k03WXt2rWIi4tDcnIywsPD8fHHH2P58uUoKipC9+7dm12fTRO1dm+t+RIZBecafX5Ib2+kTPyVjolIb598fhzz/1V83+tbABSzeXIqpy7fxOCFO+95ywFbJwu+mPEfumfSG+8Ifh8SExMxadIkTJ48GUFBQUhKSoLNZsPSpUuNjkbkFJpqmABgR0nzV+FQy/ZzGiYAqFIpB6ln1Ee7Gr1H05nr/I7djU3TD6qrq3HgwAFERkbaLY+MjERubu4916mqqkJFRYXdg6i1+tU72xWNC/nzvzROQkaJ/WS3Kq/Tc/oWVV6Hfr7Piy/g6q3bTY7h9+tHbJp+cOnSJdTW1qJrV/srPLp27Ypvvvnmnuv89a9/RadOneofNptNj6hEhjh/1zlMTam8XadxEjLKofLrRkcglRWUXTM6QovCpuknTD+5bFZEGiy7Y8aMGbh+/Xr948yZM3pEJDJE145uisZ5uPLXSmv1WLdORkcglfX1+4XREVoU/nb7wYMPPgiz2dxgVunChQsNZp/usFgs8PT0tHsQtVZ7/vs5ReN4FV3rlfZ6mCqvw6vonMczgV3wQHvXJsfw+/UjNk0/cHNzQ//+/bF9u/15G9u3b8eAAQMMSkXkXGKeaPqeX0N6e+uUhIzyl+FBP2t9i0o5SD0bEyIavUeTrRO/Y3fjLQfucueWA3//+98RFhaGTz75BMuWLUNhYSF69OjR7Pq85QC1FQPm/x/KK368qob3aWp74pbvQUHZNfT1+wXmjuqD1/93H45dbPrzyThj4dy+OHYRcSv21n/dlr5fvE/TfUpOTsZ7772Hc+fOISQkBIsXL8agQYMUrcumiYiIqOVh02QANk1EREQtD29uSURERKQiNk1ERERECrBpIiIiIlKATRMRERGRAmyaiIiIiBRg00RERESkAJsmIiIiIgXYNBEREREpwKaJiIiISAEXowO0Jndurl5RUWFwEiIiIlLqzt/t5j4khU2Tim7cuAEAsNlsBichIiIiR924cQOdOnVq9Hl+9pyK6urqUF5ejo4dO8JkMt3361RUVMBms+HMmTP8DDsFWC/HsF7KsVaOYb0cw3o5Rst6iQhu3LiBbt26oV27xs9c4kyTitq1awc/Pz/VXs/T05M/SA5gvRzDeinHWjmG9XIM6+UYrerV1AzTHTwRnIiIiEgBNk1ERERECrBpckIWiwWzZs2CxWIxOkqLwHo5hvVSjrVyDOvlGNbLMc5QL54ITkRERKQAZ5qIiIiIFGDTRERERKQAmyYiIiIiBdg0ERERESnApklj2dnZeOmll9CtWzeYTCZs2LCh2XXS0tIQGhqK9u3bw8fHB/Hx8bh8+fI9x65ZswYmkwmjRo1SN7hBtKjXypUrYTKZGjy+++47DbdEH1rtX9euXUNCQgJ8fHzg7u6OoKAgZGZmarQV+tGiXoMHD77n/jV8+HANt0R7Wu1bSUlJCAwMhNVqhc1mw5tvvsmfxUbqdfv2bcydOxcBAQFwd3dHaGgotm7dquFW6Od+6vXRRx8hKCgIVqsVgYGBSE1NbTAmIyMDwcHBsFgsCA4Oxvr161XNzaZJYzdv3kRoaCg+/PBDReNzcnIwfvx4TJo0CYWFhUhPT8e+ffswefLkBmNPnTqFadOmYeDAgWrHNoxW9fL09MS5c+fsHu7u7lpsgq60qFd1dTWee+45nDx5Ep999hmKi4uxbNky+Pr6arUZutGiXuvWrbPbr44cOQKz2Yzo6GitNkMXWtQqLS0N06dPx6xZs3D06FGsWLECa9euxYwZM7TaDN1oUa+ZM2fi448/xgcffICioiL89re/xcsvv4z8/HytNkM3jtZr6dKlmDFjBmbPno3CwkLMmTMHCQkJ2LRpU/2Y3bt3IyYmBnFxcTh48CDi4uIwZswY5OXlqRdcSDcAZP369U2OWbhwofj7+9stW7Jkifj5+dktq6mpkfDwcFm+fLlMmDBBoqKiVE5rPLXqlZKSIp06ddIgoXNRq15Lly4Vf39/qa6u1iKm01Dz5/Fuixcvlo4dO0plZaUaMZ2CWrVKSEiQoUOH2o3505/+JBEREapldQZq1cvHx0c+/PBDuzFRUVESGxurWlZnoKReYWFhMm3aNLtlU6dOlfDw8Pqvx4wZI88//7zdmGHDhsnYsWNVy8qZJiczYMAAlJWVITMzEyKC8+fP47PPPmsw1T937lx07twZkyZNMiipc1Bar8rKSvTo0QN+fn4YMWJEq/hP7X4oqdfGjRsRFhaGhIQEdO3aFSEhIZg/fz5qa2sNTG4MpfvX3VasWIGxY8eiQ4cOOiY1npJaRURE4MCBA9i7dy8A4MSJE8jMzGzxhzLvh5J6VVVVNZgRt1qtyMnJ0Tuu4Rqrxd69e3H79m0A3880RUZG2o0ZNmwYcnNz1QuiWvtFzYKCblpEJD09XTw8PMTFxUUAyMiRI+3+68/JyRFfX1+5ePGiiEibnmkSab5eu3fvllWrVklBQYFkZ2fL6NGjxWq1SklJiYbp9adWvQIDA8ViscjEiRNl//79snr1avHy8pI5c+ZomF5/atXrbnl5eQJA8vLyVE5rLDVrtWTJEnF1da0f87vf/U6j1MZRq16vvvqqBAcHS0lJidTW1kpWVpZYrVZxc3PTML3+lNRrxowZ8tBDD8n+/fulrq5O9u3bJ126dBEAUl5eLiIirq6ukpaWZrdeWlqaqvVi06QjJTtGYWGh+Pj4yHvvvScHDx6UrVu3Sp8+fWTixIkiIlJRUSE9e/aUzMzM+nXactPUXL3upba2VkJDQ2XKlCkqJzaWWvV65JFHxGazSU1NTf2yRYsWyUMPPaRVdENosX+9/vrrEhISokFaY6lVqx07dkjXrl1l2bJlcujQIVm3bp3YbDaZO3euxlugL7XqdeHCBYmKipJ27dqJ2WyW3r17yxtvvCFWq1XjLdCXknrdunVL4uPjxcXFRcxms3Tr1k3efvttASDnz58Xke+bpn/84x9263366adisVjUy6raK1GzlOwY48aNk1deecVu2RdffFHfTefn5wsAMZvN9Q+TySQmk0nMZrMcP35cwy3Qlxr1aszkyZMbHPtu6dSq16BBg+TZZ5+1G5OZmSkApKqqStXMRlJ7/7p586Z4enpKUlKS2lENp1atIiIiGpyXsmrVKrFarVJbW6tqZiOpvW99++23UlZWJnV1dfL2229LcHCw2pENpXRmTkSkurpazpw5IzU1NZKcnCwdO3as33dsNpskJibajU9MTJTu3burlpXnNDmZW7duoV07+2+L2WwGAIgIHn30URw+fBgFBQX1j5EjR2LIkCEoKCiAzWYzIrZhmqvXvYgICgoK4OPjo3k+Z6OkXuHh4Th+/Djq6urqx5SUlMDHxwdubm76hXUCjuxf//znP1FVVYVx48bpls+ZKKlVY2Pk+3/g9QnqJBzZt9zd3eHr64uamhpkZGQgKipKt5zOxtXVFX5+fjCbzVizZg1GjBhRX8ewsDBs377dbnxWVhYGDBigXgDV2i+6pxs3bkh+fn79DFFiYqLk5+fLqVOnRERk+vTpEhcXVz8+JSVFXFxcJDk5WUpLSyUnJ0eeeOIJeeqppxp9j9Z0eE6Les2ePVu2bt0qpaWlkp+fXz/F2xrOO9GiXqdPnxYPDw/5/e9/L8XFxbJ582bp0qWLzJs3T/ftU5uWP48RERESExOj27ZoTYtazZo1Szp27CirV6+WEydOSFZWlgQEBMiYMWN03z61aVGvPXv2SEZGhpSWlkp2drYMHTpUevXqJVevXtV781TnaL2Ki4tl1apVUlJSInl5eRITEyNeXl7y9ddf14/ZtWuXmM1mWbBggRw9elQWLFggLi4usmfPHtVys2nS2I4dOwRAg8eECRNE5PuG55lnnrFbZ8mSJRIcHCxWq1V8fHwkNjZWysrKGn2P1tQ0aVGvP/7xj9K9e3dxc3OTzp07S2RkpOTm5uq4VdrRav/Kzc2Vp59+WiwWi/j7+8s777xjd45TS6VVvYqLiwWAZGVl6bQl2tOiVrdv35bZs2dLQECAuLu7i81mkzfeeKNVNAFa1Gvnzp0SFBQkFotFvL29JS4uTs6ePavjVmnH0XoVFRVJ3759xWq1iqenp0RFRclXX33V4HXT09MlMDBQXF1d5dFHH5WMjAxVc5tE2ticKBEREdF94DlNRERERAqwaSIiIiJSgE0TERERkQJsmoiIiIgUYNNEREREpACbJiIiIiIF2DQRERERKcCmiYiIiJxadnY2XnrpJXTr1g0mkwkbNmxw+DVEBO+//z569+4Ni8UCm82G+fPnO/QaLg6/KxEREZGObt68idDQUMTHx2P06NH39RpTp05FVlYW3n//ffTp0wfXr1/HpUuXHHoN3hGciIiIWgyTyYT169dj1KhR9cuqq6sxc+ZMpKWl4dq1awgJCcG7776LwYMHAwCOHj2Kxx57DEeOHEFgYOB9vzcPzxERNSI1NRXe3t6oqqqyWz569GiMHz/eoFRE9FPx8fHYtWsX1qxZg0OHDiE6OhrPP/88jh07BgDYtGkT/P39sXnzZvTq1Qs9e/bE5MmTceXKFYfeh00TEVEjoqOjUVtbi40bN9Yvu3TpEjZv3oz4+HgDkxHRHaWlpVi9ejXS09MxcOBABAQEYNq0aYiIiEBKSgoA4MSJEzh16hTS09ORmpqKlStX4sCBA3jllVccei+e00RE1Air1Ypf//rXSElJQXR0NAAgLS0Nfn5+9dP+RGSsL7/8EiKC3r172y2vqqqCt7c3AKCurg5VVVVITU2tH7dixQr0798fxcXFig/ZsWkiImrCa6+9hieffBJnz56Fr68vUlJS8Jvf/AYmk8noaESE7xsis9mMAwcOwGw22z3n4eEBAPDx8YGLi4tdYxUUFAQAOH36NJsmIiI19OvXD6GhoUhNTcWwYcNw+PBhbNq0yehYRPSDfv36oba2FhcuXMDAgQPvOSY8PBw1NTUoLS1FQEAAAKCkpAQA0KNHD8XvxavniIiasXTpUixevBiRkZE4duwYtm3bZnQkojalsrISx48fB/B9k5SYmIghQ4bAy8sL3bt3x7hx47Br1y4sWrQI/fr1w6VLl/Dvf/8bffr0wYsvvoi6ujo8+eST8PDwQFJSEurq6pCQkABPT09kZWUpzsGmiYioGRUVFfDx8UFNTQ1SU1MRExNjdCSiNmXnzp0YMmRIg+UTJkzAypUrcfv2bcybNw+pqak4e/YsvL29ERYWhjlz5qBPnz4AgPLyckyZMgVZWVno0KEDXnjhBSxatAheXl6Kc7BpIiJSYPz48diyZQvKy8thsViMjkNEBuAtB4iIFDh37hxiY2PZMBG1YZxpIiJqwpUrV5CVlYXY2FgUFRX9rLsJE1HLxqvniIia8Pjjj+Pq1at499132TARtXGcaSIiIiJSgOc0ERERESnApomIiIhIATZNRERERAqwaSIiIiJSgE0TERERkQJsmoiIiIgUYNNEREREpACbJiIiIiIF2DQRERERKfD/XWLaHnBWaawAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Plot cross-section\n", + "df_red.plot.scatter(\n", + " x=\"y\", y=\"red_band_value\", title=\"Sentinel-2 red band values in y-direction\"\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "f76abba7", + "metadata": {}, + "source": [ + "The cross-section view shows most points having a Red band reflectance value of 10000,\n", + "that should correspond to white sea ice. Darker values near 0 would be water, and\n", + "intermediate values around 6000 would be thin ice.\n", + "\n", + "(Click 'Show code cell content' below if you'd like to see the histogram plot)" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "02bb268b", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [ + "hide-cell" + ] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlEAAAGxCAYAAABC0OPBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAABEcUlEQVR4nO3deVxV1f7/8feRSUA4CgSIIlJOGVqphWiphVNOlbdrhaGlv7LMgcQ0b31LG8Sh1LrezLymppXZTcvKi2KZZY6pVGaplUMWiAMCToCwfn/4YN+OoOIWBeX1fDzO49FZ+3P2XntxgrfrrL2PwxhjBAAAgPNSpbw7AAAAcDkiRAEAANhAiAIAALCBEAUAAGADIQoAAMAGQhQAAIANhCgAAAAbCFEAAAA2EKIAAABsIEQBZWzdunW6++67VadOHXl5eSkkJEQxMTFKTEy8qMc9duyYRo8erS+//LLYttmzZ8vhcGjXrl0XtQ9169bVgw8+eM66tLQ0PfPMM4qJiVFQUJD8/f3VvHlzvfnmmyooKLiofTyX0aNHy+FwnLOuXbt2ioqKugQ9Kr1L9XOuKMcFyhshCihDn332mVq1aqXs7GxNmDBBy5Yt06uvvqrWrVvr/fffv6jHPnbsmMaMGVNiiOratavWrFmjmjVrXtQ+lNbGjRv19ttvKzY2Vm+//bY+/PBDtW3bVo899pgefvjh8u4eAJSKe3l3ALiSTJgwQZGRkVq6dKnc3f/3v9d9992nCRMmlFu/rrrqKl111VXldvzTtW7dWr/++qs8PDystg4dOigvL0//+te/NGbMGIWHh5/XPvPz8+VwOFzGHQAuJmaigDJ08OBBBQUFlfiHvEqV4v+7vf/++4qJiZGvr6+qVaumTp06afPmzS41Dz74oKpVq6ZffvlFXbp0UbVq1RQeHq7ExETl5uZKknbt2mWFpDFjxsjhcMjhcFgfrZX0cUvRx1EbNmzQrbfeKh8fH1199dUaN26cCgsLXfqQnZ2t4cOHKzIyUp6enqpVq5YSEhJ09OhRW+NUo0YNlwBV5Oabb5Yk7d2796yv//LLL+VwODR37lwlJiaqVq1a8vLy0i+//CJJWr58uWJjY+Xv7y8fHx+1bt1an3/+ebH9fPbZZ7rhhhvk5eWlyMhIvfzyy+d9Ll9//bVatmwpb29v1apVS//3f/9X7CPJMWPGKDo6WgEBAfL391ezZs00c+ZMnf7973Xr1lW3bt2UnJysZs2aydvbW40aNdJbb71V7Lhr165V69atVbVqVYWFhWnUqFHKz88/Z3+nTJkih8NhjdVfjRw5Up6enjpw4IAkKSUlRXfeeadq166tqlWrql69ehowYIC1/WzO9NFuu3bt1K5dO5e20r6/PvjgA0VHR8vpdFrv1379+p2zL8DFQogCylBMTIzWrVunIUOGaN26dWf9ozZ27Fjdf//9aty4sRYsWKC5c+cqJydHt956q7Zu3epSm5+frx49eig2NlYff/yx+vXrp8mTJ2v8+PGSpJo1ayo5OVmS1L9/f61Zs0Zr1qzR//3f/521v+np6erdu7ceeOABLV68WHfccYdGjRqlefPmWTXHjh1T27ZtNWfOHA0ZMkT//e9/NXLkSM2ePVs9evQoFgQuxBdffCF3d3c1aNCgVPWjRo3Snj179MYbb+iTTz5RcHCw5s2bp44dO8rf319z5szRggULFBAQoE6dOrkEqc8//1x33nmn/Pz8NH/+fE2cOFELFizQrFmzSt3f9PR03Xffferdu7c+/vhj3XPPPXrxxRc1dOhQl7pdu3ZpwIABWrBggRYuXKiePXtq8ODBeuGFF4rt87vvvlNiYqKeeOIJffzxx2ratKn69++vr776yqrZunWrYmNjdfjwYc2ePVtvvPGGNm/erBdffPGcfX7ggQfk6emp2bNnu7QXFBRo3rx56t69u4KCgiRJv/76q2JiYjRt2jQtW7ZMzz77rNatW6dbbrmlVIGtNEr7/lqzZo3uvfdeXX311Zo/f74+++wzPfvsszp58mSZ9AOwxQAoMwcOHDC33HKLkWQkGQ8PD9OqVSuTlJRkcnJyrLo9e/YYd3d3M3jwYJfX5+TkmNDQUNOrVy+rrW/fvkaSWbBggUttly5dTMOGDa3n+/fvN5LMc889V6xfs2bNMpLMzp07rba2bdsaSWbdunUutY0bNzadOnWyniclJZkqVaqYDRs2uNT95z//MZLMkiVLrLaIiAjTt2/fMw/QWSxdutRUqVLFPPHEE+esXbFihZFk2rRp49J+9OhRExAQYLp37+7SXlBQYK6//npz8803W23R0dEmLCzMHD9+3GrLzs42AQEBpjS/GovG7+OPP3Zpf/jhh02VKlXM7t27S3xdQUGByc/PN88//7wJDAw0hYWF1raIiAhTtWpVl9ceP37cBAQEmAEDBlht9957r/H29jbp6elW28mTJ02jRo2K/ZxL0rNnT1O7dm1TUFBgtS1ZssRIMp988kmJryksLDT5+flm9+7dxc67pPfXmd4Lbdu2NW3btrWel/b99fLLLxtJ5vDhw2c9N+BSYiYKKEOBgYH6+uuvtWHDBo0bN0533nmntm/frlGjRqlJkybWxyBLly7VyZMn1adPH508edJ6VK1aVW3bti22ONzhcKh79+4ubU2bNtXu3bsvqL+hoaHWR2hn2u+nn36qqKgo3XDDDS597dSpkxwOR4kL2YsUFBS4vOb0jwmLbNq0Sb169VLLli2VlJRU6v7/7W9/c3m+evVqHTp0SH379i123M6dO2vDhg06evSojh49qg0bNqhnz56qWrWq9Xo/P79i43w2fn5+6tGjh0tbXFycCgsLXWaOvvjiC7Vv315Op1Nubm7y8PDQs88+q4MHDyojI8Pl9TfccIPq1KljPa9ataoaNGjg8jNZsWKFYmNjFRISYrW5ubnp3nvvLVW/H3roIe3du1fLly+32mbNmqXQ0FDdcccdVltGRoYeffRRhYeHy93dXR4eHoqIiJAk/fTTT6U61rmU9v110003SZJ69eqlBQsW6I8//iiT4wMXghAFXAQtWrTQyJEj9cEHH+jPP//UE088oV27dlmLy/ft2yfp1B8GDw8Pl8f7779fbM2Jj4+Pyx97SfLy8tKJEycuqJ+BgYHF2ry8vHT8+HHr+b59+/T9998X66efn5+MMWddH3PNNde4vOb5558vVrN582Z16NBB9evX15IlS+Tl5VXq/p9+tWHRuN5zzz3F+jt+/HgZY3To0CFlZmaqsLBQoaGhxfZZUtuZ/DXEnP76gwcPSpLWr1+vjh07SpJmzJihb775Rhs2bNDTTz8tSS5jLZXuZ3Lw4MEL6vsdd9yhmjVrWh9dZmZmavHixerTp4/c3NwkSYWFherYsaMWLlyoESNG6PPPP9f69eu1du3aEvttV2nfX23atNFHH31k/eOjdu3aioqK0nvvvVcm/QDs4DIW4CLz8PDQc889p8mTJ2vLli2SZK05+c9//mP9y76iCgoKkre3d4mLm4u2n8knn3xiLX6XpLCwMJftmzdvVvv27RUREaFly5bJ6XSeV99Ov59TUV/++c9/qmXLliW+JiQkxLqSLz09vdj2ktrOpCi0lfT6ojA0f/58eXh46NNPP3UJwh999FGpj3O6wMDAC+q7m5ub4uPj9dprr+nw4cN69913lZubq4ceesiq2bJli7777jvNnj1bffv2tdpLWpBekqpVq7r87IscOHDA5T1zPu+vO++8U3feeadyc3O1du1aJSUlKS4uTnXr1lVMTEyp+gWUJUIUUIbS0tJKvBdT0UcfRSGiU6dOcnd316+//lrsIym7imZwymqGoEi3bt00duxYBQYGKjIy8rxe26RJkzNuS01NVfv27VW7dm2lpKSoRo0aF9pVtW7dWtWrV9fWrVs1aNCgM9Z5enrq5ptv1sKFCzVx4kQr3OTk5OiTTz4p9fFycnK0ePFil4/03n33XVWpUkVt2rSRJOu2C0UzPNKpn9HcuXPP9/Qst912mxYvXqx9+/ZZs2EFBQXndS+yhx56SBMmTNB7772n2bNnKyYmRo0aNbK2FwXU02cGp0+fXqr9161bV99//71L2/bt27Vt2zaXYGTn/eXl5aW2bduqevXqWrp0qTZv3kyIQrkgRAFlqFOnTqpdu7a6d++uRo0aqbCwUKmpqXrllVdUrVo166qtunXr6vnnn9fTTz+t3377TZ07d1aNGjW0b98+rV+/Xr6+vhozZsx5HdvPz08RERH6+OOPFRsbq4CAAAUFBalu3boXdE4JCQn68MMP1aZNGz3xxBNq2rSpCgsLtWfPHi1btkyJiYmKjo4+r31u27ZN7du3lyS99NJL2rFjh3bs2GFtv+aaa2zd16patWr65z//qb59++rQoUO65557FBwcrP379+u7777T/v37NW3aNEnSCy+8oM6dO6tDhw5KTExUQUGBxo8fL19fXx06dKhUxwsMDNRjjz2mPXv2qEGDBlqyZIlmzJihxx57zFrX1LVrV02aNElxcXF65JFHdPDgQb388svn9bHl6Z555hktXrxYt99+u5599ln5+PjoX//613ndcqJRo0aKiYlRUlKSfv/9d7355pvFtl9zzTV66qmnZIxRQECAPvnkE6WkpJRq//Hx8XrggQc0cOBA/e1vf9Pu3bs1YcKEYj/X0r6/nn32We3du1exsbGqXbu2Dh8+rFdffVUeHh5q27Ztqc8bKFPlu64duLK8//77Ji4uztSvX99Uq1bNeHh4mDp16pj4+HizdevWYvUfffSRue2224y/v7/x8vIyERER5p577jHLly+3avr27Wt8fX2Lvfa5554rdhXZ8uXLzY033mi8vLyMJOvqqDNdnXfdddcV22/fvn1NRESES9uRI0fMM888Yxo2bGg8PT2N0+k0TZo0MU888YTLFWKlvTqvqD9nesyaNeusry+6Ou+DDz4ocfvKlStN165dTUBAgPHw8DC1atUyXbt2LVa/ePFi07RpU+Pp6Wnq1Kljxo0bV+K4lqRo/L788kvTokUL4+XlZWrWrGn+8Y9/mPz8fJfat956yzRs2NB4eXmZq6++2iQlJZmZM2eWeEVb165dSzzWX69oM8aYb775xrRs2dJ4eXmZ0NBQ8+STT5o333yzVFfnFSmq9/b2NllZWcW2b9261XTo0MH4+fmZGjVqmL///e9mz549xa4CLen9VVhYaCZMmGCuvvpqU7VqVdOiRQvzxRdflHgupXl/ffrpp+aOO+4wtWrVMp6eniY4ONh06dLFfP3116U6V+BicBhThjd5AQAAqCS4Og8AAMAGQhQAAIANhCgAAAAbCFEAAAA2EKIAAABsIEQBAADYwM02S6mwsFB//vmn/Pz8in3VBAAAqJiMMcrJyVFYWJiqVCnbuSNCVCn9+eefCg8PL+9uAAAAG37//XfVrl27TPdJiColPz8/Sad+CP7+/uXcGwAAUBrZ2dkKDw+3/o6XJUJUKRV9hOfv70+IAgDgMnMxluKwsBwAAMAGQhQAAIANhCgAAAAbWBNVhowxOnnypAoKCsq7K6gk3Nzc5O7uzm03AKAcEKLKSF5entLS0nTs2LHy7goqGR8fH9WsWVOenp7l3RUAqFQIUWWgsLBQO3fulJubm8LCwuTp6cnMAC46Y4zy8vK0f/9+7dy5U/Xr1y/zG8kBAM6MEFUG8vLyVFhYqPDwcPn4+JR3d1CJeHt7y8PDQ7t371ZeXp6qVq1a3l0CgEqDf7aWIWYBUB543wFA+eC3LwAAgA2EKAAAABtYE3UR1X3qs0t6vF3jul7S4xVp166dbrjhBk2ZMuWctXXr1lVCQoISEhIuer9O9+CDD+rw4cP66KOPLsnxyvNcAQAXHzNRAAAANhCiYMnLyyvvLgAAcNkgRFVi7dq106BBgzRs2DAFBQWpQ4cO2rp1q7p06aJq1aopJCRE8fHxOnDggPWao0ePqk+fPqpWrZpq1qypV1555byPm5OTo7i4OFWrVk1hYWH65z//6bJ90qRJatKkiXx9fRUeHq6BAwfqyJEj1vbZs2erevXqWrp0qa699lpVq1ZNnTt3VlpamlVTUFCgYcOGqXr16goMDNSIESNkjClV/6ZPn65atWqpsLDQpb1Hjx7q27evJOnXX3/VnXfeqZCQEFWrVk033XSTli9ffsZ97tq1Sw6HQ6mpqVbb4cOH5XA49OWXX1pt5xp/AEDFQYiq5ObMmSN3d3d98803GjdunNq2basbbrhB3377rZKTk7Vv3z716tXLqn/yySe1YsUKLVq0SMuWLdOXX36pjRs3ntcxJ06cqKZNm2rTpk0aNWqUnnjiCaWkpFjbq1Spotdee01btmzRnDlz9MUXX2jEiBEu+zh27JhefvllzZ07V1999ZX27Nmj4cOHW9tfeeUVvfXWW5o5c6ZWrVqlQ4cOadGiRaXq39///ncdOHBAK1assNoyMzO1dOlS9e7dW5J05MgRdenSRcuXL9fmzZvVqVMnde/eXXv27DmvsfirtLS0c44/cLmr+9Rnth9ARcPC8kquXr16mjBhgiTp2WefVbNmzTR27Fhr+1tvvaXw8HBt375dYWFhmjlzpt5++2116NBB0qkQVrt27fM6ZuvWrfXUU09Jkho0aKBvvvlGkydPtvb514XYkZGReuGFF/TYY4/p9ddft9rz8/P1xhtv6JprrpEkDRo0SM8//7y1fcqUKRo1apT+9re/SZLeeOMNLV26tFT9CwgIUOfOnfXuu+8qNjZWkvTBBx8oICDAen799dfr+uuvt17z4osvatGiRVq8eLEGDRp0XuNRZNq0aWcd/wYNGtjaLwDg4mAmqpJr0aKF9d8bN27UihUrVK1aNevRqFEjSac+vvr111+Vl5enmJgY6zUBAQFq2LDheR3zr68vev7TTz9Zz1esWKEOHTqoVq1a8vPzU58+fXTw4EEdPXrUqvHx8bEClCTVrFlTGRkZkqSsrCylpaW5HMfd3d3lXM+ld+/e+vDDD5WbmytJeuedd3TffffJzc1N0qmPNUeMGKHGjRurevXqqlatmn7++ecLmok61/gDACoWZqIqOV9fX+u/CwsL1b17d40fP75YXc2aNbVjx46L1o+i7xrcvXu3unTpokcffVQvvPCCAgICtGrVKvXv31/5+flWvYeHR7HXl3bNU2l0795dhYWF+uyzz3TTTTfp66+/1qRJk6ztTz75pJYuXaqXX35Z9erVk7e3t+65554zLs4vuqv4X/v41/ORzj3+AICKhRAFS7NmzfThhx+qbt26cncv/taoV6+ePDw8tHbtWtWpU0fSqbVC27dvV9u2bUt9nLVr1xZ7XjTj8u233+rkyZN65ZVXrOCxYMGC8zoPp9OpmjVrau3atWrTpo0k6eTJk9q4caOaNWtWqn14e3urZ8+eeuedd/TLL7+oQYMGat68ubX966+/1oMPPqi7775b0qk1Urt27Trj/q666ipJp9Y93XjjjZLksshcOvf4AwAqFj7Og+Xxxx/XoUOHdP/992v9+vX67bfftGzZMvXr108FBQWqVq2a+vfvryeffFKff/65tmzZogcffPC8v7vtm2++0YQJE7R9+3b961//0gcffKChQ4dKkq655hqdPHlS//znP/Xbb79p7ty5euONN877XIYOHapx48Zp0aJF+vnnnzVw4EAdPnz4vPbRu3dvffbZZ3rrrbf0wAMPuGyrV6+eFi5cqNTUVH333XeKi4srdjXfX3l7e6tly5YaN26ctm7dqq+++krPPPOMS825xh8AULHwz92LqLzuIG5XWFiYvvnmG40cOVKdOnVSbm6uIiIi1LlzZysoTZw4UUeOHFGPHj3k5+enxMREZWVlnddxEhMTtXHjRo0ZM0Z+fn565ZVX1KlTJ0nSDTfcoEmTJmn8+PEaNWqU2rRpo6SkJPXp0+e8j5GWlmaFvH79+unuu+8+r77efvvtCggI0LZt2xQXF+eybfLkyerXr59atWqloKAgjRw5UtnZ2Wfd31tvvaV+/fqpRYsWatiwoSZMmKCOHTta20sz/gCAisNhynIhyRUsOztbTqdTWVlZ8vf3d9l24sQJ7dy5U5GRkapatWo59RCVFe8/XE4u5FYFl9s/TFExnO3v94Xin7cAAAA2lGuIGj16tBwOh8sjNDTU2m6M0ejRoxUWFiZvb2+1a9dOP/74o8s+cnNzNXjwYAUFBcnX11c9evTQ3r17XWoyMzMVHx8vp9Mpp9Op+Pj4814fg3P7+uuvXS7PP/1RUezZs+es/byQ2xQAACqPcl8Tdd1117l8XUbRfXgkacKECZo0aZJmz56tBg0a6MUXX1SHDh20bds2+fn5STp1Y8ZPPvlE8+fPV2BgoBITE9WtWzdt3LjR2ldcXJz27t2r5ORkSdIjjzyi+Ph4ffLJJ5fwTK98LVq0KHbFWUUUFhZ21n6GhYVdus4AAC5b5R6i3N3dXWafihhjNGXKFD399NPq2bOnpFN3xw4JCdG7776rAQMGKCsrSzNnztTcuXPVvn17SdK8efMUHh6u5cuXq1OnTvrpp5+UnJystWvXKjo6WpI0Y8YMxcTEaNu2bed9o0icmbe3t+rVq1fe3Tgnd3f3y6KfAICKrdzXRO3YsUNhYWGKjIzUfffdp99++02StHPnTqWnp7tcveTl5aW2bdtq9erVkk7d4Tk/P7/YFU5RUVFWzZo1a+R0Oq0AJUktW7aU0+m0akqSm5ur7Oxsl8e5sEYf5YH3HQCUj3INUdHR0Xr77be1dOlSzZgxQ+np6WrVqpUOHjyo9PR0SVJISIjLa0JCQqxt6enp8vT0VI0aNc5aExwcXOzYwcHBVk1JkpKSrDVUTqdT4eHhZ6wtunv2sWPHSnHWQNkqet+dfhd3AMDFVa4f591xxx3Wfzdp0kQxMTG65pprNGfOHLVs2VLS/74OpIgxpljb6U6vKan+XPsZNWqUhg0bZj3Pzs4+Y5Byc3NT9erVre9u8/HxOWcfgQtljNGxY8eUkZGh6tWru6wnBABcfOW+JuqvfH191aRJE+3YsUN33XWXpFMzSX/93rCMjAxrdio0NFR5eXnKzMx0mY3KyMhQq1atrJp9+/YVO9b+/fuLzXL9lZeXl7y8vErd96J1XUVBCrhUqlevXuK6QgDAxVWhQlRubq5++ukn3XrrrYqMjFRoaKhSUlKs7xrLy8vTypUrrS9obd68uTw8PJSSkqJevXpJOvXdZFu2bNGECRMkSTExMcrKytL69et18803S5LWrVunrKwsK2iVBYfDoZo1ayo4OLjYF8sCF4uHhwczUABQTso1RA0fPlzdu3dXnTp1lJGRoRdffFHZ2dnq27evHA6HEhISNHbsWNWvX1/169fX2LFj5ePjY30Fh9PpVP/+/ZWYmKjAwEAFBARo+PDhatKkiXW13rXXXqvOnTvr4Ycf1vTp0yWdusVBt27dLsqVeW5ubvxRAwCgEijXELV3717df//9OnDggK666iq1bNlSa9euVUREhCRpxIgROn78uAYOHKjMzExFR0dr2bJl1j2ipFPfYebu7q5evXrp+PHjio2N1ezZs12CzDvvvKMhQ4ZYV/H16NFDU6dOvbQnCwAArih8d14pXczv3gGAyoLvzsOlxnfnAQAAVDCEKAAAABsIUQAAADYQogAAAGwgRAEAANhAiAIAALCBEAUAAGADIQoAAMAGQhQAAIANhCgAAAAbCFEAAAA2EKIAAABsIEQBAADYQIgCAACwgRAFAABgAyEKAADABkIUAACADYQoAAAAGwhRAAAANhCiAAAAbCBEAQAA2ECIAgAAsIEQBQAAYAMhCgAAwAZCFAAAgA2EKAAAABsIUQAAADYQogAAAGwgRAEAANhAiAIAALCBEAUAAGADIQoAAMAGQhQAAIANhCgAAAAbCFEAAAA2EKIAAABsIEQBAADYQIgCAACwgRAFAABgAyEKAADABkIUAACADYQoAAAAGwhRAAAANhCiAAAAbCBEAQAA2ECIAgAAsIEQBQAAYAMhCgAAwAZCFAAAgA2EKAAAABsIUQAAADYQogAAAGwgRAEAANhAiAIAALCBEAUAAGADIQoAAMAGQhQAAIANFSZEJSUlyeFwKCEhwWozxmj06NEKCwuTt7e32rVrpx9//NHldbm5uRo8eLCCgoLk6+urHj16aO/evS41mZmZio+Pl9PplNPpVHx8vA4fPnwJzgoAAFypKkSI2rBhg9588001bdrUpX3ChAmaNGmSpk6dqg0bNig0NFQdOnRQTk6OVZOQkKBFixZp/vz5WrVqlY4cOaJu3bqpoKDAqomLi1NqaqqSk5OVnJys1NRUxcfHX7LzAwAAV55yD1FHjhxR7969NWPGDNWoUcNqN8ZoypQpevrpp9WzZ09FRUVpzpw5OnbsmN59911JUlZWlmbOnKlXXnlF7du314033qh58+bphx9+0PLlyyVJP/30k5KTk/Xvf/9bMTExiomJ0YwZM/Tpp59q27Zt5XLOAADg8lfuIerxxx9X165d1b59e5f2nTt3Kj09XR07drTavLy81LZtW61evVqStHHjRuXn57vUhIWFKSoqyqpZs2aNnE6noqOjrZqWLVvK6XRaNSXJzc1Vdna2ywMAAKCIe3kefP78+dq0aZM2bNhQbFt6erokKSQkxKU9JCREu3fvtmo8PT1dZrCKaopen56eruDg4GL7Dw4OtmpKkpSUpDFjxpzfCQEAgEqj3Gaifv/9dw0dOlTz5s1T1apVz1jncDhcnhtjirWd7vSakurPtZ9Ro0YpKyvLevz+++9nPSYAAKhcyi1Ebdy4URkZGWrevLnc3d3l7u6ulStX6rXXXpO7u7s1A3X6bFFGRoa1LTQ0VHl5ecrMzDxrzb59+4odf//+/cVmuf7Ky8tL/v7+Lg8AAIAi5RaiYmNj9cMPPyg1NdV6tGjRQr1791ZqaqquvvpqhYaGKiUlxXpNXl6eVq5cqVatWkmSmjdvLg8PD5eatLQ0bdmyxaqJiYlRVlaW1q9fb9WsW7dOWVlZVg0AAMD5Krc1UX5+foqKinJp8/X1VWBgoNWekJCgsWPHqn79+qpfv77Gjh0rHx8fxcXFSZKcTqf69++vxMREBQYGKiAgQMOHD1eTJk2sherXXnutOnfurIcffljTp0+XJD3yyCPq1q2bGjZseAnPGAAAXEnKdWH5uYwYMULHjx/XwIEDlZmZqejoaC1btkx+fn5WzeTJk+Xu7q5evXrp+PHjio2N1ezZs+Xm5mbVvPPOOxoyZIh1FV+PHj00derUS34+AADgyuEwxpjy7sTlIDs7W06nU1lZWayPAgCb6j71me3X7hrXtQx7gsriYv79Lvf7RAEAAFyOCFEAAAA2EKIAAABsIEQBAADYQIgCAACwgRAFAABgAyEKAADABkIUAACADYQoAAAAGwhRAAAANhCiAAAAbCBEAQAA2ECIAgAAsIEQBQAAYAMhCgAAwAZCFAAAgA2EKAAAABsIUQAAADYQogAAAGwgRAEAANhAiAIAALCBEAUAAGADIQoAAMAGQhQAAIANhCgAAAAbCFEAAAA2EKIAAABsIEQBAADYQIgCAACwgRAFAABgAyEKAADABkIUAACADYQoAAAAGwhRAAAANhCiAAAAbCBEAQAA2ECIAgAAsIEQBQAAYAMhCgAAwAZCFAAAgA2EKAAAABsIUQAAADYQogAAAGwgRAEAANhAiAIAALCBEAUAAGADIQoAAMAGWyFq586dZd0PAACAy4qtEFWvXj3ddtttmjdvnk6cOFHWfQIAAKjwbIWo7777TjfeeKMSExMVGhqqAQMGaP369WXdNwAAgArLVoiKiorSpEmT9Mcff2jWrFlKT0/XLbfcouuuu06TJk3S/v37y7qfAAAAFcoFLSx3d3fX3XffrQULFmj8+PH69ddfNXz4cNWuXVt9+vRRWlpaWfUTAACgQrmgEPXtt99q4MCBqlmzpiZNmqThw4fr119/1RdffKE//vhDd955Z1n1EwAAoEJxt/OiSZMmadasWdq2bZu6dOmit99+W126dFGVKqcyWWRkpKZPn65GjRqVaWcBAAAqClshatq0aerXr58eeughhYaGllhTp04dzZw584I6BwAAUFHZClE7duw4Z42np6f69u1rZ/cAAAAVnq01UbNmzdIHH3xQrP2DDz7QnDlzLrhTAAAAFZ2tEDVu3DgFBQUVaw8ODtbYsWNLvZ9p06apadOm8vf3l7+/v2JiYvTf//7X2m6M0ejRoxUWFiZvb2+1a9dOP/74o8s+cnNzNXjwYAUFBcnX11c9evTQ3r17XWoyMzMVHx8vp9Mpp9Op+Ph4HT58+PxOGgAA4C9shajdu3crMjKyWHtERIT27NlT6v3Url1b48aN07fffqtvv/1Wt99+u+68804rKE2YMEGTJk3S1KlTtWHDBoWGhqpDhw7Kycmx9pGQkKBFixZp/vz5WrVqlY4cOaJu3bqpoKDAqomLi1NqaqqSk5OVnJys1NRUxcfH2zl1AAAASTbXRAUHB+v7779X3bp1Xdq/++47BQYGlno/3bt3d3n+0ksvadq0aVq7dq0aN26sKVOm6Omnn1bPnj0lSXPmzFFISIjeffddDRgwQFlZWZo5c6bmzp2r9u3bS5LmzZun8PBwLV++XJ06ddJPP/2k5ORkrV27VtHR0ZKkGTNmKCYmRtu2bVPDhg3tDAEAAKjkbM1E3XfffRoyZIhWrFihgoICFRQU6IsvvtDQoUN133332epIQUGB5s+fr6NHjyomJkY7d+5Uenq6OnbsaNV4eXmpbdu2Wr16tSRp48aNys/Pd6kJCwtTVFSUVbNmzRo5nU4rQElSy5Yt5XQ6rZqS5ObmKjs72+UBAABQxNZM1Isvvqjdu3crNjZW7u6ndlFYWKg+ffqc15ooSfrhhx8UExOjEydOqFq1alq0aJEaN25sBZyQkBCX+pCQEO3evVuSlJ6eLk9PT9WoUaNYTXp6ulUTHBxc7LjBwcFWTUmSkpI0ZsyY8zoXAABQedgKUZ6ennr//ff1wgsv6LvvvpO3t7eaNGmiiIiI895Xw4YNlZqaqsOHD+vDDz9U3759tXLlSmu7w+FwqTfGFGs73ek1JdWfaz+jRo3SsGHDrOfZ2dkKDw8/5/kAAIDKwVaIKtKgQQM1aNDggjrg6empevXqSZJatGihDRs26NVXX9XIkSMlnZpJqlmzplWfkZFhzU6FhoYqLy9PmZmZLrNRGRkZatWqlVWzb9++Ysfdv39/sVmuv/Ly8pKXl9cFnRsAALhy2VoTVVBQoJkzZyouLk7t27fX7bff7vK4EMYY5ebmKjIyUqGhoUpJSbG25eXlaeXKlVZAat68uTw8PFxq0tLStGXLFqsmJiZGWVlZWr9+vVWzbt06ZWVlWTUAAADny9ZM1NChQzV79mx17dpVUVFR5/x47Uz+8Y9/6I477lB4eLhycnI0f/58ffnll0pOTpbD4VBCQoLGjh2r+vXrq379+ho7dqx8fHwUFxcnSXI6nerfv78SExMVGBiogIAADR8+XE2aNLGu1rv22mvVuXNnPfzww5o+fbok6ZFHHlG3bt24Mg8AANhmK0TNnz9fCxYsUJcuXS7o4Pv27VN8fLzS0tLkdDrVtGlTJScnq0OHDpKkESNG6Pjx4xo4cKAyMzMVHR2tZcuWyc/Pz9rH5MmT5e7url69eun48eOKjY3V7Nmz5ebmZtW88847GjJkiHUVX48ePTR16tQL6jsAAKjcHMYYc74vCgsL05dffnnB66EuJ9nZ2XI6ncrKypK/v395dwcALkt1n/rM9mt3jetahj1BZXEx/37bWhOVmJioV199VTbyFwAAwBXB1sd5q1at0ooVK/Tf//5X1113nTw8PFy2L1y4sEw6BwAAUFHZClHVq1fX3XffXdZ9AQAAuGzYClGzZs0q634AAABcVmytiZKkkydPavny5Zo+fbpycnIkSX/++aeOHDlSZp0DAACoqGzNRO3evVudO3fWnj17lJubqw4dOsjPz08TJkzQiRMn9MYbb5R1PwEAACoUWzNRQ4cOVYsWLZSZmSlvb2+r/e6779bnn39eZp0DAACoqGxfnffNN9/I09PTpT0iIkJ//PFHmXQMAACgIrM1E1VYWKiCgoJi7Xv37nW5mzgAAMCVylaI6tChg6ZMmWI9dzgcOnLkiJ577rkL/ioYAACAy4Gtj/MmT56s2267TY0bN9aJEycUFxenHTt2KCgoSO+9915Z9xEAAKDCsRWiwsLClJqaqvfee0+bNm1SYWGh+vfvr969e7ssNAcAALhS2QpRkuTt7a1+/fqpX79+ZdkfAACAy4KtEPX222+fdXufPn1sdQYAAOByYStEDR061OV5fn6+jh07Jk9PT/n4+BCiAADAFc/W1XmZmZkujyNHjmjbtm265ZZbWFgOAAAqBdvfnXe6+vXra9y4ccVmqQAAAK5EZRaiJMnNzU1//vlnWe4SAACgQrK1Jmrx4sUuz40xSktL09SpU9W6desy6RgAAEBFZitE3XXXXS7PHQ6HrrrqKt1+++165ZVXyqJfAAAAFZqtEFVYWFjW/QAAALislOmaKAAAgMrC1kzUsGHDSl07adIkO4cAAACo0GyFqM2bN2vTpk06efKkGjZsKEnavn273Nzc1KxZM6vO4XCUTS8BAAAqGFshqnv37vLz89OcOXNUo0YNSaduwPnQQw/p1ltvVWJiYpl2EgAAoKKxtSbqlVdeUVJSkhWgJKlGjRp68cUXuToPAABUCrZCVHZ2tvbt21esPSMjQzk5ORfcKQAAgIrOVoi6++679dBDD+k///mP9u7dq7179+o///mP+vfvr549e5Z1HwEAACocW2ui3njjDQ0fPlwPPPCA8vPzT+3I3V39+/fXxIkTy7SDAAAAFZGtEOXj46PXX39dEydO1K+//ipjjOrVqydfX9+y7h8AAECFdEE320xLS1NaWpoaNGggX19fGWPKql8AAAAVmq0QdfDgQcXGxqpBgwbq0qWL0tLSJEn/7//9P25vAAAAKgVbIeqJJ56Qh4eH9uzZIx8fH6v93nvvVXJycpl1DgAAoKKytSZq2bJlWrp0qWrXru3SXr9+fe3evbtMOgYAAFCR2ZqJOnr0qMsMVJEDBw7Iy8vrgjsFAABQ0dkKUW3atNHbb79tPXc4HCosLNTEiRN12223lVnnAAAAKipbH+dNnDhR7dq107fffqu8vDyNGDFCP/74ow4dOqRvvvmmrPsIAABQ4diaiWrcuLG+//573XzzzerQoYOOHj2qnj17avPmzbrmmmvKuo8AAAAVznnPROXn56tjx46aPn26xowZczH6BAAAUOGd90yUh4eHtmzZIofDcTH6AwAAcFmw9XFenz59NHPmzLLuCwAAwGXD1sLyvLw8/fvf/1ZKSopatGhR7DvzJk2aVCadAwAAqKjOK0T99ttvqlu3rrZs2aJmzZpJkrZv3+5Sw8d8AACgMjivEFW/fn2lpaVpxYoVkk59zctrr72mkJCQi9I5AACAiuq81kQZY1ye//e//9XRo0fLtEMAAACXA1sLy4ucHqoAAAAqi/MKUQ6Ho9iaJ9ZAAQCAyui81kQZY/Tggw9aXzJ84sQJPfroo8Wuzlu4cGHZ9RAAAKACOq8Q1bdvX5fnDzzwQJl2BgAA4HJxXiFq1qxZF6sfAAAAl5ULWlgOAABQWRGiAAAAbCBEAQAA2ECIAgAAsIEQBQAAYAMhCgAAwAZCFAAAgA3lGqKSkpJ00003yc/PT8HBwbrrrru0bds2lxpjjEaPHq2wsDB5e3urXbt2+vHHH11qcnNzNXjwYAUFBcnX11c9evTQ3r17XWoyMzMVHx8vp9Mpp9Op+Ph4HT58+GKfIgAAuEKVa4hauXKlHn/8ca1du1YpKSk6efKkOnbsqKNHj1o1EyZM0KRJkzR16lRt2LBBoaGh6tChg3JycqyahIQELVq0SPPnz9eqVat05MgRdevWTQUFBVZNXFycUlNTlZycrOTkZKWmpio+Pv6Sni8AALhyOIwxprw7UWT//v0KDg7WypUr1aZNGxljFBYWpoSEBI0cOVLSqVmnkJAQjR8/XgMGDFBWVpauuuoqzZ07V/fee68k6c8//1R4eLiWLFmiTp066aefflLjxo21du1aRUdHS5LWrl2rmJgY/fzzz2rYsOE5+5adnS2n06msrCz5+/tfvEEAgCtY3ac+s/3aXeO6lmFPUFlczL/fFWpNVFZWliQpICBAkrRz506lp6erY8eOVo2Xl5fatm2r1atXS5I2btyo/Px8l5qwsDBFRUVZNWvWrJHT6bQClCS1bNlSTqfTqjldbm6usrOzXR4AAABFKkyIMsZo2LBhuuWWWxQVFSVJSk9PlySFhIS41IaEhFjb0tPT5enpqRo1apy1Jjg4uNgxg4ODrZrTJSUlWeunnE6nwsPDL+wEAQDAFaXChKhBgwbp+++/13vvvVdsm8PhcHlujCnWdrrTa0qqP9t+Ro0apaysLOvx+++/l+Y0AABAJVEhQtTgwYO1ePFirVixQrVr17baQ0NDJanYbFFGRoY1OxUaGqq8vDxlZmaetWbfvn3Fjrt///5is1xFvLy85O/v7/IAAAAoUq4hyhijQYMGaeHChfriiy8UGRnpsj0yMlKhoaFKSUmx2vLy8rRy5Uq1atVKktS8eXN5eHi41KSlpWnLli1WTUxMjLKysrR+/XqrZt26dcrKyrJqAAAAzod7eR788ccf17vvvquPP/5Yfn5+1oyT0+mUt7e3HA6HEhISNHbsWNWvX1/169fX2LFj5ePjo7i4OKu2f//+SkxMVGBgoAICAjR8+HA1adJE7du3lyRde+216ty5sx5++GFNnz5dkvTII4+oW7dupboyDwAA4HTlGqKmTZsmSWrXrp1L+6xZs/Tggw9KkkaMGKHjx49r4MCByszMVHR0tJYtWyY/Pz+rfvLkyXJ3d1evXr10/PhxxcbGavbs2XJzc7Nq3nnnHQ0ZMsS6iq9Hjx6aOnXqxT1BAABwxapQ94mqyLhPFABcOO4ThUut0twnCgAA4HJBiAIAALCBEAUAAGADIQoAAMAGQhQAAIANhCgAAAAbCFEAAAA2EKIAAABsIEQBAADYQIgCAACwgRAFAABgAyEKAADABkIUAACADYQoAAAAGwhRAAAANhCiAAAAbCBEAQAA2ECIAgAAsIEQBQAAYAMhCgAAwAZCFAAAgA2EKAAAABsIUQAAADYQogAAAGwgRAEAANhAiAIAALCBEAUAAGADIQoAAMAGQhQAAIANhCgAAAAbCFEAAAA2EKIAAABsIEQBAADYQIgCAACwgRAFAABgAyEKAADABkIUAACADYQoAAAAGwhRAAAANhCiAAAAbCBEAQAA2ECIAgAAsIEQBQAAYAMhCgAAwAZCFAAAgA2EKAAAABsIUQAAADYQogAAAGwgRAEAANhAiAIAALCBEAUAAGADIQoAAMAGQhQAAIANhCgAAAAbCFEAAAA2EKIAAABsIEQBAADYUK4h6quvvlL37t0VFhYmh8Ohjz76yGW7MUajR49WWFiYvL291a5dO/34448uNbm5uRo8eLCCgoLk6+urHj16aO/evS41mZmZio+Pl9PplNPpVHx8vA4fPnyRzw4AAFzJyjVEHT16VNdff72mTp1a4vYJEyZo0qRJmjp1qjZs2KDQ0FB16NBBOTk5Vk1CQoIWLVqk+fPna9WqVTpy5Ii6deumgoICqyYuLk6pqalKTk5WcnKyUlNTFR8ff9HPDwAAXLkcxhhT3p2QJIfDoUWLFumuu+6SdGoWKiwsTAkJCRo5cqSkU7NOISEhGj9+vAYMGKCsrCxdddVVmjt3ru69915J0p9//qnw8HAtWbJEnTp10k8//aTGjRtr7dq1io6OliStXbtWMTEx+vnnn9WwYcMS+5Obm6vc3FzreXZ2tsLDw5WVlSV/f/+LOBIAcOWq+9Rn5XLcXeO6lstxUf6ys7PldDovyt/vCrsmaufOnUpPT1fHjh2tNi8vL7Vt21arV6+WJG3cuFH5+fkuNWFhYYqKirJq1qxZI6fTaQUoSWrZsqWcTqdVU5KkpCTr4z+n06nw8PCyPkUAAHAZq7AhKj09XZIUEhLi0h4SEmJtS09Pl6enp2rUqHHWmuDg4GL7Dw4OtmpKMmrUKGVlZVmP33///YLOBwAAXFncy7sD5+JwOFyeG2OKtZ3u9JqS6s+1Hy8vL3l5eZ1nbwEAQGVRYUNUaGiopFMzSTVr1rTaMzIyrNmp0NBQ5eXlKTMz02U2KiMjQ61atbJq9u3bV2z/+/fvLzbLBQA4t/Ja1wRUNBX247zIyEiFhoYqJSXFasvLy9PKlSutgNS8eXN5eHi41KSlpWnLli1WTUxMjLKysrR+/XqrZt26dcrKyrJqAAAAzle5zkQdOXJEv/zyi/V8586dSk1NVUBAgOrUqaOEhASNHTtW9evXV/369TV27Fj5+PgoLi5OkuR0OtW/f38lJiYqMDBQAQEBGj58uJo0aaL27dtLkq699lp17txZDz/8sKZPny5JeuSRR9StW7czXpkHAABwLuUaor799lvddttt1vNhw4ZJkvr27avZs2drxIgROn78uAYOHKjMzExFR0dr2bJl8vPzs14zefJkubu7q1evXjp+/LhiY2M1e/Zsubm5WTXvvPOOhgwZYl3F16NHjzPemwoAAKA0Ksx9oiq6i3mfCQC4nFyOa6K4T1TlVSnvEwUAAFCREaIAAABsIEQBAADYQIgCAACwgRAFAABgAyEKAADABkIUAACADYQoAAAAGwhRAAAANhCiAAAAbCBEAQAA2ECIAgAAsIEQBQAAYAMhCgAAwAZCFAAAgA2EKAAAABsIUQAAADYQogAAAGwgRAEAANhAiAIAALCBEAUAAGADIQoAAMAGQhQAAIANhCgAAAAbCFEAAAA2EKIAAABsIEQBAADYQIgCAACwgRAFAABgAyEKAADABkIUAACADYQoAAAAGwhRAAAANhCiAAAAbCBEAQAA2OBe3h0AAFx6dZ/6rLy7AFz2mIkCAACwgRAFAABgAyEKAADABkIUAACADSwsB4DLFIvDgfLFTBQAAIANhCgAAAAb+DivAriQKfld47qWYU8AAEBpEaIAAFc8/rGKi4GP8wAAAGwgRAEAANhAiAIAALCBEAUAAGADIQoAAMAGQhQAAIANhCgAAAAbuE8UAJQjvv8OuHwxEwUAAGADM1Gwhbv/AgAqO0IULisX+tHHhQS48gqOl+NxLwQhG8DlolKFqNdff10TJ05UWlqarrvuOk2ZMkW33npreXcLl1B5BQPWvVwazJACuJQqTYh6//33lZCQoNdff12tW7fW9OnTdccdd2jr1q2qU6dOeXfPNv5o4GwIb6XH/0sAzpfDGGPKuxOXQnR0tJo1a6Zp06ZZbddee63uuusuJSUlnfP12dnZcjqdysrKkr+/f5n2jT90wP9cjh9BAmdCwC5/F/Pvd6WYicrLy9PGjRv11FNPubR37NhRq1evLvE1ubm5ys3NtZ5nZWVJOvXDKGuFucfKfJ/A5arOEx+UdxeAMnMx/mbg/BT9DC7GnFGlCFEHDhxQQUGBQkJCXNpDQkKUnp5e4muSkpI0ZsyYYu3h4eEXpY8AgCuPc0p59wBFcnJy5HQ6y3SflSJEFXE4HC7PjTHF2oqMGjVKw4YNs54XFhbq0KFDCgwMPONr7MjOzlZ4eLh+//33Mp9mvNwwFqcwDqcwDv/DWJzCOJzCOPxPacbCGKOcnByFhYWV+fErRYgKCgqSm5tbsVmnjIyMYrNTRby8vOTl5eXSVr169YvVRfn7+1f6/xmKMBanMA6nMA7/w1icwjicwjj8z7nGoqxnoIpUijuWe3p6qnnz5kpJSXFpT0lJUatWrcqpVwAA4HJWKWaiJGnYsGGKj49XixYtFBMTozfffFN79uzRo48+Wt5dAwAAl6FKE6LuvfdeHTx4UM8//7zS0tIUFRWlJUuWKCIiolz75eXlpeeee67YR4eVEWNxCuNwCuPwP4zFKYzDKYzD/5T3WFSa+0QBAACUpUqxJgoAAKCsEaIAAABsIEQBAADYQIgCAACwgRAFAABgAyGqnL3++uuKjIxU1apV1bx5c3399dfl3SXbkpKSdNNNN8nPz0/BwcG66667tG3bNpcaY4xGjx6tsLAweXt7q127dvrxxx9danJzczV48GAFBQXJ19dXPXr00N69e11qMjMzFR8fL6fTKafTqfj4eB0+fPhin6ItSUlJcjgcSkhIsNoq0zj88ccfeuCBBxQYGCgfHx/dcMMN2rhxo7W9MozFyZMn9cwzzygyMlLe3t66+uqr9fzzz6uwsNCquRLH4auvvlL37t0VFhYmh8Ohjz76yGX7pTznPXv2qHv37vL19VVQUJCGDBmivLy8i3HaJTrbWOTn52vkyJFq0qSJfH19FRYWpj59+ujPP/902ceVMBbnek/81YABA+RwODRlyhSX9go1DgblZv78+cbDw8PMmDHDbN261QwdOtT4+vqa3bt3l3fXbOnUqZOZNWuW2bJli0lNTTVdu3Y1derUMUeOHLFqxo0bZ/z8/MyHH35ofvjhB3PvvfeamjVrmuzsbKvm0UcfNbVq1TIpKSlm06ZN5rbbbjPXX3+9OXnypFXTuXNnExUVZVavXm1Wr15toqKiTLdu3S7p+ZbG+vXrTd26dU3Tpk3N0KFDrfbKMg6HDh0yERER5sEHHzTr1q0zO3fuNMuXLze//PKLVVMZxuLFF180gYGB5tNPPzU7d+40H3zwgalWrZqZMmWKVXMljsOSJUvM008/bT788EMjySxatMhl+6U655MnT5qoqChz2223mU2bNpmUlBQTFhZmBg0adNHHoMjZxuLw4cOmffv25v333zc///yzWbNmjYmOjjbNmzd32ceVMBbnek8UWbRokbn++utNWFiYmTx5ssu2ijQOhKhydPPNN5tHH33Upa1Ro0bmqaeeKqcela2MjAwjyaxcudIYY0xhYaEJDQ0148aNs2pOnDhhnE6neeONN4wxp36ZeHh4mPnz51s1f/zxh6lSpYpJTk42xhizdetWI8msXbvWqlmzZo2RZH7++edLcWqlkpOTY+rXr29SUlJM27ZtrRBVmcZh5MiR5pZbbjnj9soyFl27djX9+vVzaevZs6d54IEHjDGVYxxO/4N5Kc95yZIlpkqVKuaPP/6wat577z3j5eVlsrKyLsr5ns3ZwkOR9evXG0nWP6qvxLE40zjs3bvX1KpVy2zZssVERES4hKiKNg58nFdO8vLytHHjRnXs2NGlvWPHjlq9enU59apsZWVlSZICAgIkSTt37lR6errLOXt5ealt27bWOW/cuFH5+fkuNWFhYYqKirJq1qxZI6fTqejoaKumZcuWcjqdFWrsHn/8cXXt2lXt27d3aa9M47B48WK1aNFCf//73xUcHKwbb7xRM2bMsLZXlrG45ZZb9Pnnn2v79u2SpO+++06rVq1Sly5dJFWecfirS3nOa9asUVRUlMLCwqyaTp06KTc31+Wj5YokKytLDofD+uL7yjIWhYWFio+P15NPPqnrrruu2PaKNg6V5mtfKpoDBw6ooKBAISEhLu0hISFKT08vp16VHWOMhg0bpltuuUVRUVGSZJ1XSee8e/duq8bT01M1atQoVlP0+vT0dAUHBxc7ZnBwcIUZu/nz52vTpk3asGFDsW2VaRx+++03TZs2TcOGDdM//vEPrV+/XkOGDJGXl5f69OlTacZi5MiRysrKUqNGjeTm5qaCggK99NJLuv/++yVVrvdEkUt5zunp6cWOU6NGDXl6ela4cZGkEydO6KmnnlJcXJz8/f0lVZ6xGD9+vNzd3TVkyJASt1e0cSBElTOHw+Hy3BhTrO1yNGjQIH3//fdatWpVsW12zvn0mpLqK8rY/f777xo6dKiWLVumqlWrnrHuSh8H6dS/Klu0aKGxY8dKkm688Ub9+OOPmjZtmvr06WPVXelj8f7772vevHl69913dd111yk1NVUJCQkKCwtT3759rborfRxKcqnO+XIZl/z8fN13330qLCzU66+/fs76K2ksNm7cqFdffVWbNm06776U1zjwcV45CQoKkpubW7HEm5GRUSwdX24GDx6sxYsXa8WKFapdu7bVHhoaKklnPefQ0FDl5eUpMzPzrDX79u0rdtz9+/dXiLHbuHGjMjIy1Lx5c7m7u8vd3V0rV67Ua6+9Jnd3d6uPV/o4SFLNmjXVuHFjl7Zrr71We/bskVR53hNPPvmknnrqKd13331q0qSJ4uPj9cQTTygpKUlS5RmHv7qU5xwaGlrsOJmZmcrPz69Q45Kfn69evXpp586dSklJsWahpMoxFl9//bUyMjJUp04d63fn7t27lZiYqLp160qqeONAiConnp6eat68uVJSUlzaU1JS1KpVq3Lq1YUxxmjQoEFauHChvvjiC0VGRrpsj4yMVGhoqMs55+XlaeXKldY5N2/eXB4eHi41aWlp2rJli1UTExOjrKwsrV+/3qpZt26dsrKyKsTYxcbG6ocfflBqaqr1aNGihXr37q3U1FRdffXVlWIcJKl169bFbnOxfft2RURESKo874ljx46pShXXX7dubm7WLQ4qyzj81aU855iYGG3ZskVpaWlWzbJly+Tl5aXmzZtf1PMsraIAtWPHDi1fvlyBgYEu2yvDWMTHx+v77793+d0ZFhamJ598UkuXLpVUAceh1EvQUeaKbnEwc+ZMs3XrVpOQkGB8fX3Nrl27yrtrtjz22GPG6XSaL7/80qSlpVmPY8eOWTXjxo0zTqfTLFy40Pzwww/m/vvvL/GS5tq1a5vly5ebTZs2mdtvv73Ey1ebNm1q1qxZY9asWWOaNGlSYS5nL8lfr84zpvKMw/r16427u7t56aWXzI4dO8w777xjfHx8zLx586yayjAWffv2NbVq1bJucbBw4UITFBRkRowYYdVcieOQk5NjNm/ebDZv3mwkmUmTJpnNmzdbV5xdqnMuupw9NjbWbNq0ySxfvtzUrl37kt7i4GxjkZ+fb3r06GFq165tUlNTXX5/5ubmXlFjca73xOlOvzrPmIo1DoSocvavf/3LREREGE9PT9OsWTPrdgCXI0klPmbNmmXVFBYWmueee86EhoYaLy8v06ZNG/PDDz+47Of48eNm0KBBJiAgwHh7e5tu3bqZPXv2uNQcPHjQ9O7d2/j5+Rk/Pz/Tu3dvk5mZeQnO0p7TQ1RlGodPPvnEREVFGS8vL9OoUSPz5ptvumyvDGORnZ1thg4daurUqWOqVq1qrr76avP000+7/IG8EsdhxYoVJf5O6Nu3rzHm0p7z7t27TdeuXY23t7cJCAgwgwYNMidOnLiYp+/ibGOxc+fOM/7+XLFihbWPK2EszvWeOF1JIaoijYPDGGNKP28FAAAAiTVRAAAAthCiAAAAbCBEAQAA2ECIAgAAsIEQBQAAYAMhCgAAwAZCFAAAgA2EKAAAABsIUQAAADYQogAAAGwgRAEAANjw/wGBHIK1Pbd/8AAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "df_red.plot(\n", + " kind=\"hist\", column=\"red_band_value\", bins=30, title=\"Sentinel-2 red band values\"\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "6136bb4b", + "metadata": { + "lines_to_next_cell": 2 + }, + "source": [ + "To keep things simple, we'll label the `surface_type` of each ATL07 point\n", + "using a simple threshold:\n", + "\n", + "| Int label | Surface Type | Bin values |\n", + "|-----------|:---------------------:|:-------------------:|\n", + "| 0 | Water (dark) | `0 <= x <= 4000` |\n", + "| 1 | Thin sea ice (gray) | `4000 < x <= 8000` |\n", + "| 2 | Thick sea ice (white) | `8000 < x <= 14000` |" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "cf0052cf", + "metadata": {}, + "outputs": [], + "source": [ + "gdf[\"surface_type\"] = pd.cut(\n", + " x=df_red[\"red_band_value\"],\n", + " bins=[0, 4000, 8000, 14000],\n", + " labels=[0, 1, 2], # \"water\", \"thin_ice\", \"thick_ice\"\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "563468e0", + "metadata": {}, + "source": [ + "There are some NaN values in some rows of our geodataframe (which had no matching\n", + "Sentinel-2 pixel value) that should be dropped here." + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "id": "38c1956b", + "metadata": {}, + "outputs": [], + "source": [ + "gdf = gdf.dropna().reset_index(drop=True)" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "4815912d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
indexphoton_ratehist_wbackground_r_normheight_segment_heightheight_segment_n_pulse_seghist_mean_hhist_median_hx_atclayer_flagheight_segment_ssh_flaggeometryhist_mean_median_h_diffsurface_type
043888.6470580.2243125565264.000.60520216-67.327896-67.3059692.821440e+0700POINT (577809.769 1900049.731)-0.0219272
143898.0588240.1730395587664.500.66986816-67.209755-67.2199632.821441e+0700POINT (577808.777 1900044.162)0.0102082
243907.9411760.1904585587664.500.73693816-67.143631-67.1386642.821441e+0700POINT (577807.757 1900038.423)-0.0049672
343918.2941180.2504015587664.500.70109616-67.213936-67.2173772.821442e+0700POINT (577806.783 1900032.932)0.0034412
443927.8823530.1703195587670.500.62087216-67.279953-67.2771002.821442e+0700POINT (577805.803 1900027.401)-0.0028532
.............................................
9449138378.7500000.1485173209719.50-0.03286015-68.235619-68.2111282.827476e+0700POINT (567235.505 1840650.409)-0.0244901
9450138389.8666670.1139733322385.75-0.03839114-68.211784-68.2082212.827476e+0700POINT (567234.707 1840645.805)-0.0035631
94511383911.4615380.1306213322385.75-0.04604312-68.220390-68.2078552.827477e+0701POINT (567233.997 1840641.703)-0.0125351
94521384012.0833330.1332393322386.75-0.07062511-68.228592-68.2162402.827477e+0701POINT (567233.382 1840638.147)-0.0123521
94531384110.5000000.1368303322399.50-0.03843713-68.216888-68.2002332.827478e+0700POINT (567232.672 1840634.043)-0.0166551
\n", + "

9454 rows × 14 columns

\n", + "
" + ], + "text/plain": [ + " index photon_rate hist_w background_r_norm height_segment_height \\\n", + "0 4388 8.647058 0.224312 5565264.00 0.605202 \n", + "1 4389 8.058824 0.173039 5587664.50 0.669868 \n", + "2 4390 7.941176 0.190458 5587664.50 0.736938 \n", + "3 4391 8.294118 0.250401 5587664.50 0.701096 \n", + "4 4392 7.882353 0.170319 5587670.50 0.620872 \n", + "... ... ... ... ... ... \n", + "9449 13837 8.750000 0.148517 3209719.50 -0.032860 \n", + "9450 13838 9.866667 0.113973 3322385.75 -0.038391 \n", + "9451 13839 11.461538 0.130621 3322385.75 -0.046043 \n", + "9452 13840 12.083333 0.133239 3322386.75 -0.070625 \n", + "9453 13841 10.500000 0.136830 3322399.50 -0.038437 \n", + "\n", + " height_segment_n_pulse_seg hist_mean_h hist_median_h x_atc \\\n", + "0 16 -67.327896 -67.305969 2.821440e+07 \n", + "1 16 -67.209755 -67.219963 2.821441e+07 \n", + "2 16 -67.143631 -67.138664 2.821441e+07 \n", + "3 16 -67.213936 -67.217377 2.821442e+07 \n", + "4 16 -67.279953 -67.277100 2.821442e+07 \n", + "... ... ... ... ... \n", + "9449 15 -68.235619 -68.211128 2.827476e+07 \n", + "9450 14 -68.211784 -68.208221 2.827476e+07 \n", + "9451 12 -68.220390 -68.207855 2.827477e+07 \n", + "9452 11 -68.228592 -68.216240 2.827477e+07 \n", + "9453 13 -68.216888 -68.200233 2.827478e+07 \n", + "\n", + " layer_flag height_segment_ssh_flag geometry \\\n", + "0 0 0 POINT (577809.769 1900049.731) \n", + "1 0 0 POINT (577808.777 1900044.162) \n", + "2 0 0 POINT (577807.757 1900038.423) \n", + "3 0 0 POINT (577806.783 1900032.932) \n", + "4 0 0 POINT (577805.803 1900027.401) \n", + "... ... ... ... \n", + "9449 0 0 POINT (567235.505 1840650.409) \n", + "9450 0 0 POINT (567234.707 1840645.805) \n", + "9451 0 1 POINT (567233.997 1840641.703) \n", + "9452 0 1 POINT (567233.382 1840638.147) \n", + "9453 0 0 POINT (567232.672 1840634.043) \n", + "\n", + " hist_mean_median_h_diff surface_type \n", + "0 -0.021927 2 \n", + "1 0.010208 2 \n", + "2 -0.004967 2 \n", + "3 0.003441 2 \n", + "4 -0.002853 2 \n", + "... ... ... \n", + "9449 -0.024490 1 \n", + "9450 -0.003563 1 \n", + "9451 -0.012535 1 \n", + "9452 -0.012352 1 \n", + "9453 -0.016655 1 \n", + "\n", + "[9454 rows x 14 columns]" + ] + }, + "execution_count": 22, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "gdf" + ] + }, + { + "cell_type": "markdown", + "id": "146876cd", + "metadata": {}, + "source": [ + "### Save to GeoParquet" + ] + }, + { + "cell_type": "markdown", + "id": "c60e8b12", + "metadata": {}, + "source": [ + "Let's save the ATL07 photon data to a GeoParquet file so we don't have to run all the\n", + "pre-processing code above again." + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "fd168344", + "metadata": {}, + "outputs": [], + "source": [ + "gdf.to_parquet(path=\"ATL07_photons.gpq\", compression=\"zstd\", schema_version=\"1.1.0\")" + ] + }, + { + "cell_type": "markdown", + "id": "09f2f797", + "metadata": {}, + "source": [ + "```{note} To compress or not?\n", + "When storing your data, note that there is a tradeoff in terms of compression and read\n", + "speeds. Uncompressed data would typically be fastest to read (assuming no network\n", + "transfer) but result in large file sizes. We'll choose Zstandard (zstd) as the\n", + "compression method here as it provides a balance between fast reads (quicker than the\n", + "default 'snappy' compression codec), and good compression into a small file size.\n", + "```" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "id": "a466abcc", + "metadata": {}, + "outputs": [], + "source": [ + "# Load GeoParquet file back into geopandas.GeoDataFrame\n", + "gdf = gpd.read_parquet(path=\"ATL07_photons.gpq\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d7a04731", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "id": "04239b0e", + "metadata": {}, + "source": [ + "## Part 2: DataLoader and Model architecture\n", + "\n", + "The following parts will bring us one step closer to having a full machine learning\n", + "pipeline. We will create:\n", + "\n", + "1. A 'DataLoader', which is a fancy data container we can loop over; and\n", + "2. A neural network 'model' that will take our input ATL07 data and output photon\n", + " classifications." + ] + }, + { + "cell_type": "markdown", + "id": "b3124fe2", + "metadata": {}, + "source": [ + "### From dataframe tables to batched tensors\n", + "\n", + "Machine learning models are compute intensive, and typically run on specialized\n", + "hardware called Graphical Processing Units (GPUs) instead of ordinary CPUs. Depending\n", + "on your input data format (images, tables, audio, etc), and the machine learning\n", + "library/framework you'll use (e.g. Pytorch, Tensorflow, RAPIDS AI CuML, etc), there\n", + "will be different ways to transfer data from disk storage -> CPU -> GPU.\n", + "\n", + "For this exercise, we'll be using [PyTorch](https://pytorch.org), and do the following\n", + "data conversions:\n", + "\n", + "[`geopandas.GeoDataFrame`](https://geopandas.org/en/v1.0.0/docs/reference/api/geopandas.GeoDataFrame.html) ->\n", + "[`pandas.DataFrame`](https://pandas.pydata.org/pandas-docs/version/2.2/reference/api/pandas.DataFrame.html) ->\n", + "[`torch.Tensor`](https://pytorch.org/docs/2.4/tensors.html#torch.Tensor) ->\n", + "[torch `Dataset`](https://pytorch.org/docs/2.4/data.html#torch.utils.data.Dataset) ->\n", + "[torch `DataLoader`](https://pytorch.org/docs/2.4/data.html#torch.utils.data.DataLoader)" + ] + }, + { + "cell_type": "code", + "execution_count": 25, + "id": "cc9cf774", + "metadata": { + "lines_to_next_cell": 1 + }, + "outputs": [], + "source": [ + "# Select data variables from DataFrame that will be used for training\n", + "df = gdf[\n", + " [\n", + " # Input variables\n", + " \"photon_rate\",\n", + " \"hist_w\",\n", + " \"background_r_norm\",\n", + " \"height_segment_height\",\n", + " \"height_segment_n_pulse_seg\",\n", + " \"hist_mean_median_h_diff\",\n", + " # Output label (groundtruth)\n", + " \"surface_type\",\n", + " ]\n", + "]\n", + "tensor = torch.from_numpy( # convert pandas.DataFrame to torch.Tensor (via numpy)\n", + " df.to_numpy(dtype=\"float32\")\n", + ")\n", + "# assert tensor.shape == torch.Size([9454, 7]) # (rows, columns)\n", + "dataset = torch.utils.data.TensorDataset(tensor) # turn torch.Tensor into torch Dataset\n", + "dataloader = torch.utils.data.DataLoader( # put torch Dataset in a DataLoader\n", + " dataset=dataset,\n", + " batch_size=128, # mini-batch size\n", + " shuffle=True,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "id": "94930b5e", + "metadata": {}, + "source": [ + "This PyTorch\n", + "[`DataLoader`](https://pytorch.org/docs/2.4/data.html#torch.utils.data.DataLoader)\n", + "can be used in a for-loop to produce mini-batch tensors of shape (128, 7) later below." + ] + }, + { + "cell_type": "markdown", + "id": "79141dca", + "metadata": {}, + "source": [ + "### Choosing a Machine Learning algorithm\n", + "\n", + "Next is to pick a supervised learning 'model' for our photon classification task.\n", + "There are a variety of machine learning methods to choose with different levels of\n", + "complexity:\n", + "\n", + "- Easy - Decision trees (e.g. XGBoost, Random Forest), K-Nearest Neighbors, etc\n", + "- Medium - Basic neural networks (e.g. Multi-layer Perceptron, Convolutional neural\n", + " networks, etc).\n", + "- Hard - State-of-the-art models (e.g. Graph Neural Networks, Transformers, State\n", + " Space Models)\n", + "\n", + "Let's take the middle ground and build a multi-layer perceptron, also known as an\n", + "artificial feedforward neural network.\n", + "\n", + "```{seealso}\n", + "There are many frameworks catering to the different levels of Machine Learning models\n", + "mentioned above. Some notable ones are:\n", + "\n", + "- Easy: 'Classic' ML - [Scikit-learn](https://scikit-learn.org) (CPU-based) and\n", + " [CuML](https://docs.rapids.ai/api/cuml) (GPU-based)\n", + "- Medium: DIY Neural networks - [Pytorch](https://pytorch.org) and\n", + " [Tensorflow](https://www.tensorflow.org)\n", + "- Hard: High-level ML frameworks - [Lightning](https://lightning.ai/docs/pytorch),\n", + " [HuggingFace](https://huggingface.co/docs), etc.\n", + "\n", + "While you might think that going from easy to hard is recommended, there are some\n", + "people who actually start with a (well-documented) framework and work their way down!\n", + "Do whatever works best for you on your machine learning journey.\n", + "```\n", + "\n", + "A Pytorch model or\n", + "[`torch.nn.Module`](https://pytorch.org/docs/2.4/generated/torch.nn.Module.html) is\n", + "constructed as a Python class with an `__init__` method for the neural network layers,\n", + "and a `forward` method for the forward pass (how the data passes through the layers).\n", + "\n", + "This multi-layer perceptron below will have:\n", + "- An input layer with 6 nodes, corresponding to the 6 input data variables\n", + "- Two hidden layers, 50 nodes each\n", + "- Output layer with 3 nodes, for 3 surface types (open water, thin ice,\n", + " thick/snow-covered ice)" + ] + }, + { + "cell_type": "code", + "execution_count": 26, + "id": "b7a29219", + "metadata": {}, + "outputs": [], + "source": [ + "class PhotonClassificationModel(torch.nn.Module):\n", + " def __init__(self):\n", + " super().__init__()\n", + " self.linear1 = torch.nn.Linear(in_features=6, out_features=50)\n", + " self.linear2 = torch.nn.Linear(in_features=50, out_features=50)\n", + " self.linear3 = torch.nn.Linear(in_features=50, out_features=3)\n", + "\n", + " def forward(self, x: torch.Tensor) -> torch.Tensor:\n", + " x1 = self.linear1(x)\n", + " x2 = self.linear2(x1)\n", + " x3 = self.linear3(x2)\n", + " return x3" + ] + }, + { + "cell_type": "code", + "execution_count": 27, + "id": "b835cac4", + "metadata": { + "lines_to_next_cell": 2 + }, + "outputs": [ + { + "data": { + "text/plain": [ + "PhotonClassificationModel(\n", + " (linear1): Linear(in_features=6, out_features=50, bias=True)\n", + " (linear2): Linear(in_features=50, out_features=50, bias=True)\n", + " (linear3): Linear(in_features=50, out_features=3, bias=True)\n", + ")" + ] + }, + "execution_count": 27, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "model = PhotonClassificationModel()\n", + "# model = model.to(device=\"cuda\") # uncomment this line if running on GPU\n", + "model" + ] + }, + { + "cell_type": "markdown", + "id": "6e5faf45", + "metadata": {}, + "source": [ + "## Part 3: Training the neural network\n", + "\n", + "Now is the time to train the ML model! We'll need to:\n", + "1. Choose a [loss function](https://pytorch.org/docs/2.4/nn.html#loss-functions) and\n", + " [optimizer](https://pytorch.org/docs/2.4/optim.html)\n", + "2. Configure training hyperparameters such as the learning rate (`lr`) and number of\n", + " epochs (`max_epochs`) or iterations over the entire training dataset.\n", + "3. Construct the main training loop to:\n", + " - get a mini-batch from the DataLoader\n", + " - pass the mini-batch data into the model to get a prediction\n", + " - minimize the error (or loss) between the prediction and groundtruth\n", + "\n", + "Let's see how this is done!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3bfafafb", + "metadata": {}, + "outputs": [], + "source": [ + "# Setup loss function and optimizer\n", + "loss_bce = torch.nn.BCEWithLogitsLoss() # binary cross entropy loss\n", + "optimizer = torch.optim.Adam(params=model.parameters(), lr=0.001)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f7ff6874", + "metadata": { + "lines_to_next_cell": 2 + }, + "outputs": [], + "source": [ + "# Main training loop\n", + "max_epochs: int = 3\n", + "size = len(dataloader.dataset)\n", + "for epoch in tqdm.tqdm(iterable=range(max_epochs)):\n", + " for i, batch in enumerate(dataloader):\n", + " minibatch: torch.Tensor = batch[0]\n", + " # assert minibatch.shape == (128, 7)\n", + " assert minibatch.device == torch.device(\"cpu\") # Data is on CPU\n", + "\n", + " # Uncomment two lines below if GPU is available\n", + " # minibatch = minibatch.to(device=\"cuda\") # Move data to GPU\n", + " # assert minibatch.device == torch.device(\"cuda:0\") # Data is on GPU now!\n", + "\n", + " # Split data into input (x) and target (y)\n", + " x = minibatch[:, :6] # Input is in first 6 columns\n", + " y = minibatch[:, 6] # Output (groundtruth) is in 7th column\n", + " y_target = torch.nn.functional.one_hot(y.to(dtype=torch.int64), 3) # 3 classes\n", + "\n", + " # Pass data into neural network model\n", + " prediction = model(x=x)\n", + "\n", + " # Compute prediction error\n", + " loss = loss_bce(input=prediction, target=y_target.to(dtype=torch.float32))\n", + "\n", + " # Backpropagation (to minimize loss)\n", + " loss.backward()\n", + " optimizer.step()\n", + " optimizer.zero_grad()\n", + "\n", + " # Report metrics\n", + " current = (i + 1) * len(x)\n", + " print(f\"loss: {loss:>7f} [{current:>5d}/{size:>5d}]\")" + ] + }, + { + "cell_type": "markdown", + "id": "66358b37", + "metadata": { + "lines_to_next_cell": 2 + }, + "source": [ + "Did the model learn something? A good sign to check is if the loss value is\n", + "decreasing, which means the error between the predicted and groundtruth value is\n", + "getting smaller." + ] + }, + { + "cell_type": "markdown", + "id": "709f7a1b", + "metadata": { + "lines_to_next_cell": 2 + }, + "source": [ + "## References\n", + "- Koo, Y., Xie, H., Kurtz, N. T., Ackley, S. F., & Wang, W. (2023).\n", + " Sea ice surface type classification of ICESat-2 ATL07 data by using data-driven\n", + " machine learning model: Ross Sea, Antarctic as an example. Remote Sensing of\n", + " Environment, 296, 113726. https://doi.org/10.1016/j.rse.2023.113726" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5014ed28", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "jupytext": { + "formats": "py:percent,ipynb" + }, + "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.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/book/tutorials/photon_classifier.py b/book/tutorials/machine-learning/photon_classifier.py similarity index 99% rename from book/tutorials/photon_classifier.py rename to book/tutorials/machine-learning/photon_classifier.py index 21b0162..05f07bd 100644 --- a/book/tutorials/photon_classifier.py +++ b/book/tutorials/machine-learning/photon_classifier.py @@ -1,7 +1,7 @@ # --- # jupyter: # jupytext: -# formats: py:percent +# formats: py:percent,ipynb # text_representation: # extension: .py # format_name: percent