mirror of
https://github.com/heyarne/earth-observation-for-journalism.git
synced 2026-05-06 19:13:40 +02:00
340 lines
33 KiB
Text
340 lines
33 KiB
Text
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Calculating Spectral Indices Over a Long Time Span"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"from pathlib import Path\n",
|
|
"from multiprocessing import Pool\n",
|
|
"import geopandas as gpd\n",
|
|
"from sentinel_helpers import scihub_band_paths\n",
|
|
"\n",
|
|
"base_path = Path('resources/tempelhofer_feld/')\n",
|
|
"input_files = list(base_path.glob('*.zip'))\n",
|
|
"area_of_interest = gpd.read_file(base_path / 'tempelhofer_feld.geojson')\n",
|
|
"ndvi_path = Path('output') / 'ndvi'"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import rasterio as r\n",
|
|
"import rasterio.mask\n",
|
|
"import rasterio.plot as rplt\n",
|
|
"import numpy as np\n",
|
|
"from sentinel_helpers import scihub_cloud_mask\n",
|
|
"from zipfile import BadZipFile"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# uncomment this cell to remove previous calculations\n",
|
|
"# ! rm -rf {ndvi_path / '*.tif'}"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 5,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def normalize(v):\n",
|
|
" return (np.clip(v, 0, 10_000) / 10_000).astype('f4') # ← four-byte-float / float32\n",
|
|
"\n",
|
|
"def calculate_ndvi(product_path):\n",
|
|
" global area_of_interest\n",
|
|
" \n",
|
|
" try:\n",
|
|
" b04_path, b08_path = scihub_band_paths(product_path, ['B04', 'B08'], '10m')\n",
|
|
" cloud_mask = scihub_cloud_mask(product_path)\n",
|
|
" except BadZipFile:\n",
|
|
" # This is due to a problem of Scihub\n",
|
|
" print(f'{product_path}: Problem reading product, skipping it')\n",
|
|
" return\n",
|
|
"\n",
|
|
" with r.open(b04_path, 'r') as b04, r.open(b08_path, 'r') as b08:\n",
|
|
" # we want to only write the bare minimum data necessary to disk\n",
|
|
" out_meta = b04.meta.copy()\n",
|
|
"\n",
|
|
" # we reproject the geojson file we fetched above and convert it so that rasterio\n",
|
|
" # can use it as a mask; we subtract the cloud mask to avoid irregularities\n",
|
|
" reprojected = area_of_interest.to_crs(out_meta['crs']).iloc[0].geometry\n",
|
|
" mask = reprojected.difference(cloud_mask)\n",
|
|
" \n",
|
|
" if mask.area / reprojected.area < 0.25:\n",
|
|
" print(f'{product_path}: Area is covered by clouds, skipping it')\n",
|
|
" return\n",
|
|
" \n",
|
|
" miny, minx, maxy, maxx = mask.bounds\n",
|
|
" \n",
|
|
" # update the dimensions and save as geotiff, not jp2\n",
|
|
" out_meta.update({\n",
|
|
" 'width': maxx - minx,\n",
|
|
" 'height': maxy - miny,\n",
|
|
" 'driver': 'GTiff',\n",
|
|
" 'dtype': 'float32'\n",
|
|
" }) \n",
|
|
" \n",
|
|
" out_name = Path(b04_path).name.replace('B04', 'NDVI').replace('.jp2', '.tif')\n",
|
|
"\n",
|
|
" ndvi_path.mkdir(exist_ok=True, parents=True)\n",
|
|
" with r.open(ndvi_path / out_name, 'w+', **out_meta) as dst:\n",
|
|
" # we take only the part out of our source raster that we actually need\n",
|
|
" # crop=True clips off the borders\n",
|
|
" b04, transform_b04 = rasterio.mask.mask(b04, shapes=[mask], crop=True, nodata=-999)\n",
|
|
" b08, _ = rasterio.mask.mask(b08, shapes=[mask], crop=True, nodata=-999) # we ignore the returned transform because it's identical to the previous one\n",
|
|
"\n",
|
|
" b04 = normalize(b04) # <- f4 = float32\n",
|
|
" b08 = normalize(b08)\n",
|
|
"\n",
|
|
" # uncomment the following line to see if your masked shape looks correct\n",
|
|
" #rplt.show(b04, transform=transform_b04)\n",
|
|
" #rplt.show(b08, transform=transform_b04)\n",
|
|
"\n",
|
|
" # we want to be able to ignore divide by zero errors so the formula is nicer to write\n",
|
|
" np.seterr(divide='ignore', invalid='ignore')\n",
|
|
" ndvi = (b08 - b04) / (b08 + b04)\n",
|
|
"\n",
|
|
" # uncomment the following line to see if we calculated the index correctly\n",
|
|
" # rplt.show(ndvi, transform=transform_b04)\n",
|
|
"\n",
|
|
" dst.write(ndvi)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 6,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"application/vnd.jupyter.widget-view+json": {
|
|
"model_id": "5a94deac202b4ff3ad9ecb350267dfc3",
|
|
"version_major": 2,
|
|
"version_minor": 0
|
|
},
|
|
"text/plain": [
|
|
"HBox(children=(HTML(value=''), FloatProgress(value=0.0, max=36.0), HTML(value='')))"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
},
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"resources/tempelhofer_feld/S2A_MSIL2A_20190626T102031_N0212_R065_T33UUU_20190626T125319.zip: Area is covered by clouds, skipping it\n",
|
|
"resources/tempelhofer_feld/S2A_MSIL2A_20190603T101031_N0212_R022_T33UUU_20190603T114652.zip: Problem reading product, skipping it\n",
|
|
"resources/tempelhofer_feld/S2A_MSIL2A_20190404T101031_N0211_R022_T32UQD_20190404T174806.zip: Problem reading product, skipping it\n",
|
|
"resources/tempelhofer_feld/S2A_MSIL2A_20190216T102111_N0211_R065_T33UUU_20190216T130428.zip: Problem reading product, skipping it\n",
|
|
"resources/tempelhofer_feld/S2B_MSIL2A_20190419T101029_N0211_R022_T33UUU_20190419T132322.zip: Problem reading product, skipping it\n",
|
|
"resources/tempelhofer_feld/S2A_MSIL2A_20190407T102021_N0211_R065_T33UUU_20190407T134109.zip: Problem reading product, skipping it\n",
|
|
"resources/tempelhofer_feld/S2B_MSIL2A_20190512T102029_N0212_R065_T33UUU_20190512T134103.zip: Problem reading product, skipping it\n",
|
|
"resources/tempelhofer_feld/S2A_MSIL2A_20190613T101031_N0212_R022_T33UUU_20190614T125329.zip: Problem reading product, skipping it\n",
|
|
"resources/tempelhofer_feld/S2A_MSIL2A_20190424T101031_N0211_R022_T32UQD_20190424T162325.zip: Problem reading product, skipping it\n",
|
|
"resources/tempelhofer_feld/S2A_MSIL2A_20190822T101031_N0213_R022_T32UQD_20190822T143621.zip: Problem reading product, skipping it\n",
|
|
"resources/tempelhofer_feld/S2A_MSIL2A_20190623T101031_N0212_R022_T33UUU_20190623T132509.zip: Problem reading product, skipping it\n",
|
|
"resources/tempelhofer_feld/S2B_MSIL2A_20190320T101029_N0211_R022_T33UUU_20190320T195148.zip: Area is covered by clouds, skipping it\n",
|
|
"\n",
|
|
"CPU times: user 226 ms, sys: 151 ms, total: 377 ms\n",
|
|
"Wall time: 24.2 s\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"%%time\n",
|
|
"from tqdm.notebook import tqdm\n",
|
|
"\n",
|
|
"# we parallelize the NDVI calculation using a Python threadpool\n",
|
|
"with Pool() as pool:\n",
|
|
" for _ in tqdm(pool.imap_unordered(calculate_ndvi, input_files), total=len(input_files)):\n",
|
|
" # this loop is only here so we can get the progress bar\n",
|
|
" pass"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"How many files could we process?"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 7,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"25\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"! ls -l {ndvi_path} | wc -l"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Plot the NDVI\n",
|
|
"\n",
|
|
"See [02c Multithreading.ipynb](./02c Multithreading.ipynb) for a performance comparison of single vs multi-threaded iteration to average the ndvi."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 8,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"from numpy import ma\n",
|
|
"output_files = list(ndvi_path.glob('*.tif'))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 9,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def average(file_path):\n",
|
|
" with r.open(file_path) as src:\n",
|
|
" # masked=True makes sure to respect the nodata value and reads the data\n",
|
|
" # as a numpy MaskedArray, which lets us use `numpy.ma`'s methods\n",
|
|
" data = src.read(1, masked=True)\n",
|
|
" return file_path, ma.average(data)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 10,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"application/vnd.jupyter.widget-view+json": {
|
|
"model_id": "d29e20cc277a4897b076e0a275161fa5",
|
|
"version_major": 2,
|
|
"version_minor": 0
|
|
},
|
|
"text/plain": [
|
|
"HBox(children=(HTML(value=''), FloatProgress(value=0.0, max=24.0), HTML(value='')))"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
},
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"\n",
|
|
"CPU times: user 59.6 ms, sys: 104 ms, total: 164 ms\n",
|
|
"Wall time: 703 ms\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"%%time\n",
|
|
"with Pool() as pool:\n",
|
|
" averages = [avg for avg in tqdm(pool.map(average, output_files), total=len(output_files))]"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"We use the pandas library to plot the year-long NDVI development:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 14,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"<AxesSubplot:>"
|
|
]
|
|
},
|
|
"execution_count": 14,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
},
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYMAAAD2CAYAAAA0/OvUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABDFElEQVR4nO2debzcVN24n+ltoaUtFNrKjlREr1EBW0Xkhw61sorsPSyWfQkgCLGCWsUFXhFEGhGBphQBQYEjgvIKKKIQRUGliGJD4QUKUlq2lgJl6Xbn98dJLukwS3InmSTn5vl8+pl77yQz36c5M9+cvVKr1SgpKSkpGdwMyTqAkpKSkpLsKZNBSUlJSUmZDEpKSkpKymRQUlJSUkKZDEpKSkpK0CAZuK67e9YxJI2OTgG6uunqFaCzn65ucb0KnwyATbMOIAV0dArQ1U1XrwCd/XR1i+WlQzJ4MOsAUkBHpwBd3XT1CtDZT1e3WF46JIO9sg4gBXR0CtDVTVevAJ39dHWL5aVDMrgq6wBSQEenAF3ddPUK0NlPV7dYXkPTiqKLzACsrINIGB2dAnR109UroNB+c+fOHQIcBowC1lqDZ5111hFz586VmQSWDhVgOfBRYlyzSrk2UUlJie7MnTv388A/Jk2a9FjWsXSDuXPnvg/42KRJk34W9ZzCNxO5rntd1jEkjY5OAbq66eoVoIHfqGaJQAO3dzBp0qTHVq5c+cU455Q1g5LEcGzv3cAhqCrqq8ArDf69alrGmsyCLBmUzJ0798RJkybNzjqObhLXufB9Bq7rXletVqdlHUeSFNjpNGB6u4Mc21tOKDnQOGm8I4n4j8+blrE8jeA7ocDXLBI6++nqtnLlymOByMmgrBmUJIZjez8FdgU+BmxQ92/9Bn9r9m94i7dZARxrWsbPU5Eo0ZI81AwqlUoNmFmr1ab7v38ZGFWr1b5dqVS+DZwAvAiMBB4GvlGr1bxKpXI1cF+tVnNCr7U/cGKtVtu7Uqksr9Vqo+rfbzDWDC6vVqsnZx1HkhTYaSzwomkZzwPPNzogiptje+vQPJEcC/zMb5I637SMXNzNFPiaRUJnvy66rQAOrFQq36vVai81eN6u1Wo/AKhUKocAf6xUKh8Grge+CjihYw/1/96UlStXHk6MmkHhkwFwZtYBpEBRncYBjQp5mLZupmWsRN0hvVj/nGN7N6DGT58HbO3Y3hdMy1g9gFiTpqjXLCo6+3XLbTXqy9kCvt7qwFqtdmOlUvkscDjwY+DqSqWyaa1WW1ypVNYDPoOqSTRl2LBhv4wTnA7JYDrwnayDSJiiOo0FnmhzTEdupmWscGxvGvAU8DVgC8f2DslBP0JRr1lUtPFzbO+HwA7B78PW3ejd8x/0nu7wZR8yLeOMCMddCvy7Uql8P8KxDwK9tVptTaVSuRkQwMXAvsDdtVrttVYnr1q1ajfgkgjvAySYDIQQe6IC7QHmSCnPb3Lcx4D7gUOklDfFObcJOrYdF9VpLLCkzTEdu5mW0QfMcGzvaeAywHVsbx/TMhZ3+todUNRrFhVt/VavGvJCt96rVqu9WqlUfgp8EXizzeGV0M/XAxeivicPBX7a7r16enr+Hie2RJKBEKIHlfF2AxYC/xBC3Cql9BocdwHwu7jntmAy8H+dW+SKwjk5tjcUGEP7ZqLE3EzLcBzbWwjcCNzv2N5epmVELTdJU7hrFhNt/Orv4F3XPbFarXazc/mHqLv+dstFfAR4wP/5L8CmlUple2BnVEJoSV9f3/vjBJXUpLMdgcellE9KKVcCNwD7NTjuNOCXwAsDOLcZ8wcYc54potNG/mO7mkGibqZl3AZUgXWAvzq2t2uSrx+DIl6zOOjs11W3Wq22FJDAcc2OqVQqBwG743cS19SwTwlcA9xeq9Xeavc+lUrluThxJdVMtDnwTOj3hcDHwwcIITYHDgA+jRp6GPlc//wTgRMBent7H5kyZcoVqDvRXVzXfQvYCbgFOAU4G5hTrVaPdl332mq1eoTrulcDxwPnopoWDkA1V01AVdeWAb3A3ahOm4uAC6vV6snBOOTQo43qwDwGuAOYCARNFJuisv5eqMw/o1qtWg1e43JUx9V0VBV8MqpQFtIJxn8XoGdo32uu636rgdMIYAEwzXXdBUk69U7k+PkPjv/E0GFrHlq9qufOm+Xf/36g2HGXLlynwGknYLnrunvl/TrFdAqXvW8BU4rqNGLEiK1d1/0iDT5PqL6QyWk7DRkypMffcGbTGTNm3HbBBRecMXLkyPX8Y18ZPXr0OZVKZdpGG230HuAPP/zhD+/ZYYcd3gw+T1/96ldXnH/++dsffPDBv3Jd95DgOlUqlYrrut+rL3vDhw9/pd6pWq2Gv2vXplardfxvqmJO6Pcjpk6dekndMb+YOnXqTv7PV0+dOvXgqOe2+nfPPfcckoRDnv4V0WnWzHmfnDVzXm3WzHm7ZeU2a+a8DWfNnHePH8fXZs2cVymvWelXq9V44IEHTtTVrdm/++67b3ac45NqJloIbBn6fQtgUd0xHwVuEEI8BRwMXCaE2D/iua1YEDfYAlBEp7H+Y7tmotTcTMt4GdgDdWd4HjDL78voBkW8ZnHQ2U9LtyFDhrTrv1uLpD4o/wC2FUJMAJ5FdW4cHj5ASjkh+FkIcTXwGynlr4QQQ9ud24adgFi95gWgiE7j/Md2BTBVN3/o6RHA06ihp8NQE9XSpojXLA46+2np1tfX9544xydSM5BSrgZORY0SekT9Sc4TQpwkhDhpIOfGePtbBhh2nimiU9SaQepupmX0mZYxA/gRcKRje5un/Z4U85rFoeh+lRbPFd2tIT09PQ/FOT6xKrSU8nbg9rq/zWpy7NHtzo3BKag7QJ3IhZNjezuhRi8calrGP9scPhY13f6NNsd10+1HqPHcx6I60NIkF9csRYrut3zu3Lnva7KMddHd3sHcuXPft3r16p3jnKPDDOSzsw4gBTJ3cmzvQ6gEvSFq6Ga7ZDAOeCnCWkFdczMt4wnH9u4CTnBs77yUl87O/JqlTNH9rgcOmzt37mTqdjobOXLk03Pnzj0xm7BSoQIs7+vrOyjOSTokgznA0VkHkTCZOjm2NwG4EzWcrgf4QITTosw+hu67OcAvUB3LA619RkHHchim0H6TJk3qAxru+uW67tX+EFOt8IeVHh31+HIJ65K1cGxvE+BeVI3gU6gv0zWmZVTbnHcvsMK0jCnpRxkdx/aGoeax/M20jDiTGUtKBhU6bHt5bdYxJE1WTo7tjUF15G8C7G1axjzUxJ3eCKePI0LNoNtupmWsQk1W2ifNjmQdy2EYnf10dYvrVdYMSgBwbG89VNPQx4DPmpZxl//3L6MWyBprWsbSFue/CNxkWkbu1rx3bO89qNVUv2VaxjlZx1NSkkd0qBlcnXUMSdNtJ38zmZuATwCHB4nA5xH/sWntwLG9Iai1idpOcsniepmW8SQq0R3v2F5PGu+hYzkMo7Ofrm5xvQqfDFDrbuhG15z8L8drUGu/mKZl1G+IESzi1aqpaAyqLEXpQM7qes1GzXTfM6XX17EchtHZLxduju19wrG9bRN8yVheOiSDtMePZ0FXnBzbq6A2vzgUOMu0jDkNDnsKWEnrZBBMOIsy/T2r63UraitOM6XX17EchtHZL3M3/7N4a8KxxHotHZLBZVkHkALdcjoHOBm4wLSMCxsd4I/Nf4zWw0uDpSii1AwyuV5+R/JPgM86trdlu+MHgI7lMIzOfnlwex/qc7RFgq8Zy0uHZHBA1gGkQOpOju1ZwDeAK2g/+/IRotUMoiSDLK/XFagJOWmsVaRjOQyjs18e3ILZwpsl+JqxvHRIBvdnHUAKpOrk2N5RwExUp/HJEWYNzwfe49jeuk2ej9NMlNn1Mi1jAW93JCc94VLHchhGZ788uAXJYFO/ySgJYnnpkAwmtD+kcKTm5NjefsCVwO+BaRGXaJiPKivvbfJ8nGairK/XbFRVPOmO5Ky90kZnvzy4BclgOGpARhLE8tIhGbTbVLqIpOLkbwl5I2rJ8QNNy1gR8dRgRFGzfoOxwGrg1QivlfX1+l/gOZLvSM7aK2109svUzbG9DQEDeNj/06Ytjv2mf0MXhVheOiSDZVkHkALLkn5Bx/Y+ihqt8DhqUtnyGKc/6j826zcYByyJ0NwEGV8vvyP5SmDvhDuSlyX4WnlkWdYBpMiyjN//E/7jTf5jw34Dx/ZGo/r5don4usviBKFDMoiyVELRSNTJsb1e1D6sS4DdW80kboRpGa8D/20RV9RF6mjxGt1kDqojuemG5AMgD15porNf1m47A2tQN2vQvGawG2qzptsivm4sLx2Swd1ZB5ACiTk5trcVqn9gDbCbaRlxthQN02qNorFE6zyGHFwv0zKeQq3BlGRHcuZeKaOzX9ZuOwMPoWrt0HxE0T7AK8BfIr5uLC8dkkGcLTKLQiJOju2NR42eGY2qETze5pRWzAd6/aUn6om0SJ1PXq7XbGBzYO+EXi8vXmmhs19mbv7NyMeBv/pNt6/RoGbgf+72Bn7rN3VGIZaXDsngoqwDSIGOnRzbWx/4LbAVqo/g3x2+5CPASNQXaD1xmonycr1+AywGktrUJC9eaaGzX5Zu2wPrAX/1f19M45rBJGBjVLmNSiwvHZJBw5mzBacjJ8f2RqDaH7cDDjItI2q1shUN1yjyx0THaSbKxfUKdSTv5TeldUouvFJEZ78s3YIhpcFndBGN+wz2AfpQN3hRieVV+GRQrVZzt2Ryp3Ti5Fc7b0BtTHOkaRl3JBRWswXrRqM6tSLVDHJ2vRLrSM6ZV+Lo7Jex287AQtMynvF/X0zzZHC/aRlRb7piexU+Gbiue13WMSTNQJ38dsUrgX2BL5iWcX2CYT2P6ryqn2sQZymKXF0v0zKeRt1pddyRnCevNNDZL2O3nXm7iQhUzWCz8Cxkx/Y2AyYSr4kotlfhk0G1Wp2WdQxJMxAnv/DMBI4EzjYt4/IkY/LnEDRaoyiYfRzpjiWH12s2qo32s528SA69EkVnv6zcHNvbAtWnF04Gi4ERwPqhvwWDHGIlg7hehU8GOt6xDNDpG8DpwA+B7yYa0Ns0Gl5a2JqBz29Qd2MddSTn0CtRdPbL0C3oL6ivGcDancj7oOb5/CfOi5c1Aw2I6+TY3imo5aivAaZHnAk8EOajFtLaIPS3OIvU5e56mZaxmrc7kt890NfJm1fS6OyXodvOqCUjHgr9bbH/uCmAY3vDUZPNfhP3cz0YawZ21jEkTRwnx/YOA36MGj10vGkZfakF1rgTOc4idXm9Xlf6jwPuSM6pV2Lo7Jeh287A3+vmDdTXDHZFDT2N1UQE8b0KnwyA87IOIAUiOTm2tzfwU8AFDvHvctOk0X7IY4Ea0ddByd318juS7wCO66AjOXdeCaOzX9fdHNtbD/gIazcRQV3NANWX9SZwzwDeJpaXDsngmKwDSIG2To7t7QL8EvgXsJ9pGW+lHhUsAFbxzmSwNOJS2JDf6+Wg7sb2GeD5efVKCp39snD7KDCUumRgWsZrwOu8PaJoH+Au0zIGsrJqLC8dkkFS4+jzREsnx/a2R1Ubnwb2Mi0jytLRHeNXZx/nnc1EUWcfQ36v1+3Aswy8IzmvXkmhs18Wbu/xH70GzwUTzwxgawbQROQTy0uHZDAx6wBSoKmTY3vvRS2y9ipqvaEXuxaVYj5rzzWIM/sYcnq9Qh3Jezq2t/UAXiKXXgmis18WbkFfW6PPbzDxLKilRl2ltJ5YXolt/SeE2BO4GOgB5kgpz697fj/gXNSU6tXAGVLKe/3nnkIt0LQGWC2l/GiMt17c/pDC0dDJsb3NUSuQ9qBWIP1vV6NSPAJ8zrG9YX5NYRxq2FtU8ny9rkQN0T3ef4xDnr2SQGe/LNzGASuBRvuKLEI1I+0D/NO0jGcH+B6xvBKpGQgheoBLgb1QVZvDhBBG3WF/ALaXUu6A2pB8Tt3zk6WUO8RMBIMGx/bGolYgHQvsaVrGo21OSYv5qJuIbfzf4yxSl2v85HoHcKxje8OyjqdEa8YBLzUZLroYNRltZwZeK4hNUjWDHYHHpZRPAgghbgD2I9QeJqUMZ8CRqBEoSdB0i7gCs5aTY3ujUIViG1QimJtJVIrw8NL5xG8myvv1clDDdPcBbolxXt69OkVnvyzcxtH8c7MIWMf/eaD9BRDTK6lksDnwTOj3hag1utdCCHEA8D3gXaw9/b8G3CmEqAGOlHJ2g3NPxO/c6+3tfWTKlClXoDaOHu+67o7ATqgP7ynA2cCcarV6tOu611ar1SNc170aVf0/F7gMOAC4H7Vp9JuooZG9qA0hDkct/3phtVo92XXd66rV6rTQo40atnUM6k5yImsPCXsQVUu6CphRrVatBq9xOXAmMB34OTAZ9eW6ltNbbwy9bZ3ho+9a+VbPlhtvudzd/+CP35Ol08ZbDv/F88+MBvjAH//4p0Ng3Igx49/8DEAbpxGo0Ui9rutumdfr9P4dXjz00YfGPzt8vVU/Bm6J6LQT4Lmu+708OsUse2Gn8OfpAOAazZzmVKvVo4GDfbeuOQ0ZMnarnmG19VzXHV/vNH6z5fu/uGgUQ4b0vbl17zLPdd1vxXVyXfda4Pv1TtVqNfw9vTa1Wq3jf1MVc0K/HzF16tRLWhz/qalTp94V+n0z//FdU6dO/dfUqVM/FfW977nnnrOScMjTv8Bp1sx5PbNmzrtp1sx5tVkz5x2ZdVzBv1kz5y2cNXPeNbNmztvSj+0Ena7XrJnzvj1r5ry+WTPnba2TVyf/8uo3a+a8XWfNnDehaG6zZs6bP2vmvBubPPdp/3N1VTe9khpNtBAIby6+BW/PpHsHUso/AdsIIcb5vy/yH19AZbkdY7z3VbGjzT9X+WOMHeAgwDIt46cZxxQmWKMo1lIUPkW4XleiaqsnxDinCF4tcWxvgxZLcuTOz7G9qcAfgWs7fKks3FoNyZ4PrAA6XXU4lldSyeAfwLZCiAlCiHWAQ3l7c2cAhBDvFUJU/J8notrElgghRgohRvt/HwnsTrwFmWYkIZAzZgDno5ZH+B/TMn6YbTjvIBheGmspCp/cXy9/bfnbideRnHuvVvjLn98JPNxks59c+Tm2VwWuQ41C/H+O7dUPWIlDV90c2+sBNqLJTZS/T/n6pmXc2eFbxfJKJBlIKVcDp6LGvz+i/iTnCSFOEkKc5B92EPAfIcRDqJFHh0gpa6it3O4VQvwL+Dtwm5Qy8m4+1WrVSsIhT8x/cPxi4CxUu+U3Mw6nEfNRm9p82P89zoYbRbleDrAJ8LkoBxfIqxnHoGrkI4HZ4fX0IV9+ju19EPgV8CSqb3IV8Wpxa5GB24aoTZWafm5My1jZ6ZvE9arUamktcNkdgk6XrONICsf2TkCtsX89MC3lhecGhGN7U4C7gJuBA4FNTMt4Psq5Rble/hpFC4B5pmXs2e74ong1wrG9McBjwKPAjcAlwDGmZVwdHJMXP38PgPtQc20+YVrG047t3Qh8Bth8IMuydNvNsb1e1E3z503L+Hla7xPXq/AzkPNQQJPCsb2DgVmo0QdH5TER+ATDS/+f/7g06olFuV7+jOQ5wO6O7U1od3xRvJrwHVT/z2mo2uiTqN3y+smDn790+u3ABsDe/gKDAFegml0OHMjrZuAWa0OogTIYl7DWYtMNx/Z2Qw0fu2/b7V96pW5Z27yxCDVzcmMgVqwFu16RO5IL5tWPY3sfAr4AOKZlPOTfgDyL+nLtJ2s/x/bWRQ0u+QBwoGkZD4We/iMqgQ1oXakM3LqSDMrNbQqIY3sfRxV0D9jn05/+1GEZh9QSf9ZkUDuINfu4SNfLtIyFqMl+bTuSi+QV4PcLXILa2/rs0FNLqUsGWfr5ndtXo8bZH2Naxl3h5/0EdgVQdWzvfXFfv6wZKAqfDPzJJoXFvzO7A3gONbt4WUGcgr0NYhXogriFcVA1oH1bHVRAL4CpqM1Tvm5aRjipvyMZZOx3AWqE4ldNy2h2t3s1as2z2B3JGbh1q2YQy6vwyQA167CQ+G3Rd6JmN+5mWsZz/lNFcBpQzYBiuIX5LWp2vdnmuEJ5ObY3EjWD9p+ou+ow70gGZOTn2N7pwJdRIxC/3+w4/7NzK3C036QUh267jQPeNC3jjZTfJ5aXDslgetYBDATH9jZBrUA6HNjDtIwFoaeL4DTQZFAEt378TXvmALs5tveeFocWygv4Gmpy6GkNNiZaAoxwbG9E6G9d9/MnldmoJtTTI+wBPBv1RbtfzLfqtlurdYmSJJaXDskgtaFZaeEP5fstahz73qZl1E+yK4JTkAziFuoiuNVzJWrp9VZNEIXxcmxvG9Rd43WmZfylwSHB6LANQ3/rqp9je59CTSr7K2oIZpSd9H6P2vApbkdyt69dt5JBLC8dksHkrAOIg7/36W9QS30fYFrG/Q0OK4LT46iNOea3O7COIrithb+e/G9QHcnrNDmsSF42ai39s5o8HySDsaG/dc3Pn1T2a9QIof2ibvnodyTPAab4CS8q3b52cVf6HSixvHRIBnG/jDLD/yK5CfgEcLhpGb9vcmjunfwZklvzzvbmduTerQkOarXdZh3JhfBybG8v1Kzqc0zLaLb5SZAMwv0GXfHzN3D6LfAWakvXuM2QV6E2yTo+xjndvnbdqhnE8tIhGYzJOoAohIbH7QWcZFrGTS0OH9ONmDrFtIw3BjAxbkwasXSB36F2dGvWkTyme6EMDL9j9WLUTOOLWxzaKBmMSSmsfvxJZXegJpXtZVrGU3Ffw6/F3QYcE2NdqTFx36dDupUMxsQ5WIdkMKL9IdkSGs99GPAV0zLa3U3n3qkDCukW6kj+TJMmiCJ4nQFsi+qMbbX2TaNkkKpf3aSyg+omlcVlNmo4cKR1pejitfMT1Bi6kwxieemQDBa0PyRzvoPafOL7pmU0HR4XoghOA6XIbj9BNUE06kjOtZff/HI28GvTMn7X5vCgaSacDFLz82vNV6HauI9t0Xwald+iltWP2pHczWsX/J92IxnE8tIhGeyUdQCtcGzvDPwdiICvRjwt104dUli3UEfyMQ06kvPudSFqZ8MvRTj2ddRKoOFkkKbfBaha89dMy+h0b4KgFnclal2prSOc0s1r15UJZz6xvHRIBnH2qe0qju0diRq58UtUP0HUJWJz65QARXcLOpLrx7Ln1ssfpnkYqmb6ZLvj/XJaP/EsFb+6SWUXJPjSP0GtK3VchGO7ee26mQxieemQDE7JOoBGOLa3L6pA3kX0cdIBuXRKiKK73Ykay17fkZxLL38p7ktQnd/nxzh1KWsPLU3cz1+lN86kssiYlvFfVGf0cf7/QSu6ee26mQxieemQDM5uf0h3cWxvV0ACD6DmEqyI+RK5c0qQQruFOpKnOLb33tBTefUyge2A6TGXP6ivGSTqF5pUdh/xb5aicgWwKfDZNsd189p1MxnE8tIhGczJOoAwju1NQq2R8gTwWdMylg/gZXLllDA6uDXqSM6dl2N744BzUUs8/zLm6fXJIDG/0KSyp4B9o04qGwC3AYtpv3hdN6/dQLaKHSixvAq/01me8Hcw+jNqrf9d/A7HEg1xbO8W1OY+WySxRWEaOLbnoNrMtzctY17Mc68CppiW0Wg/5E5i2hxVGxiG2qnsqSRfv8H7nYvaC3iC33SUKY7t2cBxpmWsn3Us9RS+ZuC6bsejD5LA30T8TtQd426dJIK8OKWBRm4OMB7YH/Ln5ddQTwAuiZsIfNaqGSTh59jexqh2/A1Ra3I91elrRuBK1H7DxzY7oMvXrlsTzmJ7lTWDBHBsbzyqRrAJUDUt418Zh1SSMv7Y+CeBJ0zLmJJ1PGH82O4FtgHeZ1rGKwN4ja8D/wMMH0CfV/h1RqJGXn0e2N3/894JzCWIE8NvgQ8CW6fUNxEnljuAsaZl7JhlHI3QoWZwdZbv79je+qhJLlsB+ySRCLJ2ShNd3EK7a33asb1tc+Y1DbX+1VcHkgh81pqFHMfPsb1hju3t5djedcDzwM+AD6P2T9ium4nAZzZque49Gz3Z5WvXzZrB1XGObzfkqgjEWZAqURzbG47qCNsOtbrivQm9dGZOXUAnt58A30Y1x+TCy785+T7wN+CaDl4qnAwW08bPX3Ll46hEJFBNaC+jEsHPgHsHsI5VUvwvKimdiOpUrqeb124cb+8SmDaxvApfM0CNlug6/tjlG4EqcJRpGbcn+PKZOHUJbdz8VT9vBY5Zs6ZyXtbx+HwTNSnu1A6/fOvXJ2p33Q5BdQwfB9yNahraxLQM07SMP2WYCDAtYxVquYvP+h3Y9XSzTHatZkBMLx1qBpd1+w39Ntk5qOWMTzUtI+nNMbru1EV0c5sNHPjfx8Y8xaeTf3HH9t6PGpDw4wjHfgA4HbjStIwHOnzr+vWJ2l23d/uPW5mW8WKH750GwXIwx6D6QsJ0pUz6LQmj6F4yiOWlQ83ggG6+mV8dvgg4CvimaRmXpvA2XXXqMrq5/R54qm9N5YyUXv8HwCV+809T/HJ5MWpY84wE3re+ZtDuuo1ELf/QrS+6WJiW8QTwB+B4/2YuTLfKZDCju1v/R7G8dEgGjXYKS5Ovo5YCvph33mEkRbeduolWbkFH8qqVPds6tve+JF/bX2QtmD27RZvD9wd2Q92gJHFnXp8M2l23kcAbSS4pkQKzUTWY3er+3q0y2c3ZxxDTS4dkMKFbb+TY3smodrhrgS+lWPC75pQBOrr9BGrNlrbuBBM1Rh5gy2YH+Vup2sDDwOUJvfdrqDkzwd1su+s2ErXaaZ75FWqr1vqlrbtVJrudDGJ56ZAM0prKvhaO7R2GWlnxVtQMwjQ7xLrilBHauZmW8dywdfr+Bhztb9LSMX778vHAXP9PTZMBai/jdwOnmZaxOon3b7Byabvrlvtk4M8UvwbY17G9TUJPdatMdjsZxPJKrANZCLEnqumkB5gjpTy/7vn9UHfVfcBq4Awp5b1Rzm3Dss6jb41je3sDPwX+BBzij05Ik2Upv36WLMs6gDRYb/TKm19ZMuIHqHbaGxJ4yYNRXx5HooZDNkwGflPSV4AbTctwE3jfMOFksKzNsblPBj5XoJbMPpq3V3Fd1qX37nYyWBbn4ERqBkKIHtRd816AARwmhDDqDvsDsL2UcgfU1PA5Mc5tRW9n0bfGsb1dUJvY/xu1qNZbab6fT6pOGaOl2yZbLV+O2lmq2R7JcTkFeAy19/JimtcMZqJusL6c0PuGWcLbyaDddStEMjAt4zHgHuAkx/aC/9NulckgGSxteVRyxPJKqploR+BxKeWTUsqVqDujtTb/kFIul1IGbezByINI57bh7s5Cb45je9ujdrZ6BrVB96tpvVcdqTnlAC3dKhX+iLrr3NUfDjpg/D2WPwE4fnPkMzRIBo7t7YaqiXzXtIyFnbxnE8I1g3bXrRDJwOd7qKVjHnNs79zXXxv2ty697zjg5aSa8iIQ67OWVDLYHFVgAxb6f1sLIcQBQoj5qGrvsXHObcHh8UKNhr9W/e9QHWm7mZbxQhrv04RUnHKCrm6HoyY2rabzjuRN/cf/+I/vSAb+tps/Qi2VflGH79eMcDJod90KkwxMy7gTeD+qQ/kbzz6x/p8d2zu6wZDTpOnmhDOI+VlLqs+g0uBv7xhpI6W8BbhFCPEpVP/BZ6KeK4Q4EX8UQG9v7yNTpky5AhgDLHBdd0fUfp+3oKrXZwNzqtXq0a7rXlutVo/w1+k43n/fy1B3VPejetzfRLWv9QJ3v/XGUBPGHDxkSG2Dvr4hH+md+OJ5wDTXda+rVqvTXNe1gfNQE1juACaiqvKgPsgPopq9rgJmVKtVK3Ru8Hg5cCYwHfg5ajPw+Wk5oQrGRcCF1Wr15AbxpO00AtWM8rrrultq5rQTcEPvxBdPn//guF8PGVI71bG9b/ROfPGKgTgNHbbRv1ev6gF4zXXd62D8M5VK7QDH9iq9E1+cCZy3/kajf/Lq0uG9G238hv2uzV+vuq6buFPP0LEr+9ZUNnVddyiwLaiVMBs5DekZu/mQIbUXXNc9JOfX6RbglN6JnA2smP/g+P83dFjfrStXDLlq2DprfjD74nlT37/DS7unUfZ6ho59X63GKtd1j0rDibrPE3By/XWqVqvhG++1SGTVUiHEJ4BvSyn38H//GoCU8nstzlkAfAxVyGKdG8Z13cur1erJnRm8jWN7G6E6it8NTE5gJmdsknbKE7q6BV5+082dwOGmZVw/kNdybE+gljr5sGkZ/3Fsz0L1DYw1LWOpY3ubAo8CfzYto90uXgPGsb2zgXOAdXsnvnhxq+vm2N6zwB2mZeRijaY43HOPe/mj/xx/L6pDeQvgZuAsf6JaYji290/gGdMy9k3ydZsR97OWVM3gH8C2QogJwLPAodRVUYQQ7wWekFLWhBATgXVQHVTL2p3bioQTwSjgduC9qD6CricCSNYpb+jqFvL6A2ppaxMYUDIARvuPr/mPQX/Alqimm/OBdVGTH9Mk6OjcMMJ1K0wzUT277lo9eddd+zcsmo5atmIfx/YuRvXHDHTl13rGAv9M6LXaEvezlkgbmZRyNXAqqo39EfUnOU8IcZIQ4iT/sIOA/wghHkKNHjpESllrdm7U91bV6M7xx4ffDHwUNXw0s47OpJzyiK5ugZff4TsbqPo73w2E+mQQVO23dGxvZ9Rw05mmZfzfQOONSP8s5FbXzV8KYxQFTQaha/eGaRnnolorfo5qovk/x/ZO8hem7JRxdGe7SyD+Z63c3AZwbK8HNYrpYOBo0zI6Wfq3ZJDj7+i1EPiRaRnTB3B+0DyzjmkZq/yVNheibpqOBTYGege4v3acOHZH3aTtYlrGX1octy7wFvAN0zK+m2ZM3cTfLc4GPonqzP/SQPdi8GeJv47aY+KC5KJMjsLPQO70TtO/q5mFSgRfykMi0PXuGfR1C3uZlvE8qlPvaH8mcVxGAytCkxufQ41Smo7qXD0z7UTgE6lmgGoigoLXDOoxLWMuaon6g1GOdzq295sB1vi6PeEs9met8MmgWq1O6/Alvofqbf+uaRl2AiF1TAJOuUVXtwZes1HDMg8awMuN5u0mIvytGhehRrX8iWRmOEehPxm0uW6FTgat3EzLqJmW8UvgA6hlPz4FPOzY3o8c2xvb7LwGdD0ZxP2sFT4Z+EO4BoRje2eipvJfjhqWlQs6cco7uro18Pojag5A/aJoUVgrGfg8g5pp/MUurgwaJIOxba5boZNBlDJpWsYK0zIuRA0uuRL4Aqo/4XTH9oZFeJssagaxPmuFTwaosbyxcWzveNT2gDegNqjJU+dJXnbNSgNd3dbyCnUkf8rfdCYOjZLBRcBJSeyxHYNXUQloI1pft0InA2KUSdMyXjAt4yRgB+AB4IfAfxzb+5zf5NyMricDYn7WdEgGx8Q9wbG9gwAHtZH9UVluydeE2E4FQle3Rl5XA6uIPyP5HcnAtIxbTMu4YmChDQz/c/EyKhm0um5BMuhGP0YaxC6TpmU8DOwB7IOaJHsr8HvH9rZrckoWySCWlw7J4I44Bzu29xnUsLH7gYP8ZW3zRiyngqGr2zu8/CVMbgGOitmR3KhmkBXBYnWtrlvRawYDKpN+f8JtwIdR241OBP7p2J7j2N676g4fh6plLesk0JjE8tIhGUyMeqBjezui1iOZD+xjWsYbaQXVIZGdCoiubs28HNSX6cExXms0qokmDwTrE7W6bkVPBh2VSdMyVpmW8SNUf8IlqOG/jzu295XQTcA4YKk/GKBbxPLSIRksbn9IPy8Afwb2MC3j5ZTiSYI4TkVDV7dmXvcAjxOvIzlPNYMgGbS6bkVPBomUSdMylpqWcQbwIcBFzRT3HNs7GBhP9/eHjuWlQzKIjGkZT5mWsZdpGc9lHUvJ4CDUkfxJx/ai7tORt2TQbghl0ZNBopiW8ahpGZ9D7bW8HPgFcCDdTwax0CEZbNr+kMKho1OArm6tvK4hYkdyaGmHPCWDjWjtV/RkkEqZNC3jLuAjqHWqlqD2qO4msbx0SAYPZh1ACujoFKCrW1MvvyP5ZlRH8og2rzMStax7npLB+qtWDvl3i2OCZJDXPrh2pFYmTctYY1rGbNRmOqem9T5NiOWlQzLYK+sAUkBHpwBd3dp5OcCGtO9Irl+kLmuWAPStqezf4piRwJs5HKIdldTLpGkZfRn8/8Ty0iEZXJV1ACmgo1OArm7tvO5Bjdnfuc1xeUsGSwFeXDSy3dDSojYRweAtk2uhQzKYkXUAKaCjU4Cubi29/BnurwHt5hvkMhmsv9GKM1ocU/RkMCjLZD2FTwbVatXKOoak0dEpQFe3iF5vUdBksGjB+te2OKbQyWCQl8l+Cp8MdFwSWUenAF3dInoVNhmM22x5q7vMwm5sA4O+TPZT+GSg45LIOjoF6OoW0auwyeClRaN+0eKYotcMBnOZ7KfwyUDHrK6jU4CubgnWDNb3H/OSDJYBtQ3Gvnlsi2MKnQwGeZnsp/DJQMesrqNTgK5uutYMgpVLX1ky4k8tDit0MhjkZbKfwicD13UvzzqGpNHRKUBXt4hecZJBnr5cl44YtXKPFs8XOhkM8jLZz9C0AukiZ2YdQAro6BSgq1sUr6jJYHnOJnAtfXP5sFarqBY6GTC4y2Q/ha8ZoDYJ1w0dnQJ0dYviFTUZ5KKJKMTSocP6Wi2yV/RkMJjLZD86JIOfZx1ACujoFKCrWxSvwiaD1auHNFyH39//dxjFTgaDuUz2o0MymJx1ACmgo1OArm5RvAqbDCrNl7Eu+paXMLjLZD86JIP5WQeQAjo6BejqFsWrqMlgSa1WWc+xvZ4GzxV9+WoY3GWyHx2SwZisA0iBMVkHkCJjsg4gJcZEOOYtYLi/Z0Ez8pgMlvqPYxo8p0MyGJN1ACkxJs7BOiSDduvDFxEdnQJ0dYvi9Zb/uE6LY/KcDDZq8JwOyWAwl8l+dEgGC7IOIAV0dArQ1S2KV5AMWjUV5TkZNOo30CEZDOYy2Y8OyWCnrANIAR2dAnR1i+IVNRm0GtOfBbrXDAZzmewnsUlnQog9gYuBHmCOlPL8uuc/D3zF/3U5cLKU8l/+c0+h7obWAKullB+N8da3dBh6HtHRKUBXtyheLZOBY3tD/efyWjPQNRkM5jLZTyI1AyFED3Apaps1AzhMCFE/SWUBUJVSbgecC8yue36ylHKHmIkA4JSBxJxzdHQK0NUtile7mkGu1iUKscR/1DUZDOYy2U9SNYMdgcellE8CCCFuAPYDvOAAKeVfQ8ffD2yR0HufndDr5AkdnQJ0dYviVdRksMx/1DUZDOYy2U9SyWBz4JnQ7wuBj7c4/jggvKdqDbhTCFEDHCllfa0BIcSJwIkAvb29j0yZMuUK1NCpL7muexaqfewWVDY8G5hTrVaPdl332mq1eoTrulcDx6NqJZcBB6CS0gTgTVSB7wXuBg4HLgIurFarJ7uue121Wp0WerSB84BjfI+JwGI/1E2BB1G1pKuAGdVq1WrwGpej1g6ZjpopOBk1LlhXpxGo2uElruserJnTTsAurus+0crpXVuMOPSFhaPYZKvXjnRd98Z6pzHjRj2z7KUR9AztW+m67uU5cLoFOKV3Imc/9tDYvr6+IRvVO20wdtTUV5aMYJsPLZnmuq5bgOv0js8T8CiwTYHLXsPvCGCNuiRvO1Wr1fD39FpUarVas+ciI4SYCuwhpTze//0IYEcp5WkNjp2M+o/eRUq5xP/bZlLKRUKIdwG/B06TUrZaMrekpHA4tjcZ+COwq2kZboPndwLuA/Y2LaPVBvRdx7G9J4D7TMuYVvf3bwHfBoaaltFwyYqSYpDUaKKFwJah37cAFtUfJITYDpgD7BckAgAp5SL/8QVUltsx6hv7GVArdHQK0NUtoldRm4kYtu6aDWg+tHRFkRPBIC+T/STVTPQPYFshxATgWeBQVDWqHyHEVsDNwBFSysdCfx8JDJFSvub/vDtwTtQ3rlarRyQQf67Q0SlAV7eIXoVNBqtW9MyleZ9BkfsLBnuZ7CeRmoGUcjVwKvA74BH1JzlPCHGSEOIk/7Bvou4sLhNCPCSEeMD/+8bAvUKIfwF/B26TUv426nv77XxaoaNTgK5uEb0KmwzWG7VyWzRNBoO8TPaT2DwDKeXtwO11f5sV+vl4VEdG/XlPAtt38NbveE0N0NEpQFe3KF6FTQZvLB92B6rGX0/hkwGDu0z2o8MM5HOzDiAFdHQK0NUtildhk8HoMSu2BzZ0bK/+O0OHZDCYy2Q/OiSDy7IOIAV0dArQ1S2KV5RksBpYkUhECfLWG8PuAirABnVP6ZAMBnOZ7EeHZHBA1gGkgI5OAbq6RfGKkgxeMy2j8/HeCTNqg5Xb+D/WjyjSIRkM5jLZjw7J4P6sA0gBHZ0CdHWL4hXc8bdMBsmEkyyrVw150P+xvhNZh2QwmMtkPzokgwlZB5ACOjoF6OrW1su0jD5gJQVMBsNHrgqWndAxGQzaMhlGh2TwZtYBpICOTgG6ukX1arX15frkNBkMHdrXbLG6kRR7/2MoyySgRzJYlnUAKbAs6wBSZFnWAaTEsojHtUoGua0ZDB3Wt9D/UceawbKsA0iJZXEO1iEZ9GYdQAro6BSgq1tUr0ImgxGjVm3u/9ifDBzb6wHWpfjJYLCXSUCPZHB31gGkgI5OAbq6RfUqZDIYMoQ/oHZgC48m0mH5aijLJKBHMji8/SGFQ0enAF3donoVMhmg/JaydjORLslgsJdJQI9kcFHWAaSAjk4BurpF9WqYDBzbq5DvZHARasczHZPBYC+TgB7J4MKsA0gBHZ0CdHWL6tWsZjActX94XpPBhehbMxjsZRLQIBlUq9WTs44haXR0CtDVLYZXs2QQrEv0ajIRJYvvV58MRvmPhU4GZZlUFD4ZuK57XdYxJI2OTgG6usXwapcMclkz8P20rBmUZVJR+GRQrVantT+qWOjoFKCrWwyvQiYD328psFFo5VItkkFZJhWFTwY6ZnUdnQJ0dRtENYMhqJnSoEkyKMukovDJQMesrqNTgK5ug6hmAG83FWmRDMoyqSh8MnBd1846hqTR0SlAV7cYXoVMBr5f/fpEWiSDskwqCp8MgPOyDiAFdHQK0NUtqlchkwHKT8uaAWWZBPRIBsdkHUAK6OgUoKtbVK+3gOH+JLMweU8Gx9A4GawyLWNVNiElxmAvk4AeyeCOrANIAR2dAnR1i+oVLCu8Tt3fg2SQ1+Wg76BxMih6rQDKMgnokQwmZh1ACujoFKCrW1SvZltfjgbeNC1jdXIhJcpE4GX/52CxOl2SwWAvk4AeyWBx1gGkgI5OAbq6RfVqlQzy2kQEsNi0jJWomotuNYPBXiYBPZJBSUmRKGoyCAgvVqdLMihBj2SwadYBpICOTgG6ukX1KmoyCPzCS1LokgwGe5kE9EgGD2YdQAro6BSgq1tUr6Img8CvPhnktcM7DoO9TAJ6JIO9sg4gBXR0CtDVLapXUZNB4KdjzWCwl0lAj2RwVdYBpICOTgG6ukX1KmoyCPyWot9oosFeJgE9ksGMrANIAR2dAnR1i+pV1GQQ+AUrl1bQJxkM9jIJwNCk3lUIsSdwMWq3pjlSyvPrnv888BX/1+XAyVLKf0U5txXVatVKIPxcoaNTgK5uMbwKmQxCfktQn9PRqM1tCp8MyjKpSKRmIIToAS5FtVEZwGFCCKPusAVAVUq5HXAuMDvGuU3RcflZHZ0CdHWLuYQ1hJKBvz/AKHKcDEJ+wSzk8SiHwieDskwqkqoZ7Ag8LqV8EkAIcQOwH+AFB0gp/xo6/n5gi6jntkLH5Wd1dArQ1S2G1wtAH7A9cKP/t2D7yNwmg5BfkAyCz2/hk0FZJhVJJYPNgWdCvy8EPt7i+ON4e92MSOcKIU4ETgTo7e19ZMqUKVcAY4AzgenATsAtwCnA2cCcarV6tOu611ar1SNc170aOB5VK7kMOACVlCag1otZBvQCdwOHAxcBF1ar1ZNd172uWq1OCz3aqBUBj/E9JvL2bL9NUUO69kJ14MyoVqtWg9e4PBT7z4HJwHyNnUagaoeXAgdq5rQTsCvwaDun3olc9sS8DR9dvbLnmHvudv9dGcIbm01YZ+iiBRuw7vDVI1zX/VaOnMJl79FqtbrNZhNeOWXRgg0Yv9nyc19cNIoNx7/5Gdd1byjQdWr0eXqyWq2+p8Blr+F3BFABVoedqtVq+Lt2LSq1Wq3Zc5ERQkwF9pBSHu//fgSwo5TytAbHTkb9R+8ipVwS59ySEh1wbO8w/A+3aRn3OLbXCzwCfN60jJ9nG11rHNv7IPAfVOfkecDRpmVck21UJUmQ1GiihcCWod+3ABbVHySE2A6YA+wnpVwS59xm+NlTK3R0CtDVLabXr1FNQkf6vwcrlr6aaFAJEvILmomCz2zhm4nKMqlIqpnoH8C2QogJwLPAoahqVD9CiK2Am4EjpJSPxTm3DWd2EnhO0dEpQFe3yF6mZbzh2N5NwMGO7Z1K/vcygLf9gmSwlf9Y+GRAWSaBhGoGUsrVwKnA71DVXSmlnCeEOEkIcZJ/2DdRk1UuE0I8JIR4oNW5Md5+ehIOOUNHpwBd3eJ6XYtKAvtRjGQwHcC0jBWoBKBNzYCyTAIJzjOQUt4O3F73t1mhn49HdWREOjcGuW5jHSA6OgXo6hbXy0UNnDgidG6ek0HYbyl6JYOyTKLHDOTJWQeQAjo6BejqFsvLtIw+4DpgD+C9/p/znAzCfkuBDf2fdUgGZZlEj2QwP+sAUkBHpwBd3QbidS3qM3ic/3uek0HYb2noZx2SQVkm0SMZjMk6gBQYk3UAKTIm6wBSYkzcE0zLeAR4ADWCrg94I+GYkmRM6GfdksGYrANIiTFxDtYhGYzIOoAU0NEpQFe3gXpd6z8uNy2j80k/6RH2WxL6WYdkUJZJ9EgGC7IOIAV0dArQ1W2gXjegZonmuYkI1vYLagZrgJUZxJI0ZZlEj2SwU9YBpICOTgG6ug3Iy7SMF1CT0BYmG07ihP2CZPB6zmszUSnLJAkOLc2QW7IOIAV0dArQ1a0Tr6OAdZIKJCXCfkEy0GHLSyjLJKBHzeCUrANIAR2dAnR1G7CXaRmvm5bxcpLBpEDYr79mkEUgKVCWSfRIBmdnHUAK6OgUoKubrl4BYT/dkoGu1y6Wlw7JYE7WAaSAjk4Burrp6hUQ9tMtGeh67WJ5JbKEdUlJyeDBsb3NUItK3mlaxh5Zx1OSDIWvGfibOGiFjk4Burrp6hVQ5xf0b2hRM9D12sX1KmsGJSUlsXFs7w3gl6ZlHJF1LCXJoEPN4OqsY0gaHZ0CdHXT1Suggd/fgIczCCVxdL12cb10mGfQcFnsgqOjU4Cubrp6BazlZ1qGTit96nrtYnkVvmaA2uhZN3R0CtDVTVevAJ39dHWL5aVDMrgs6wBSQEenAF3ddPUK0NlPV7dYXjokgwOyDiAFdHQK0NVNV68Anf10dYvlpUMyuD/rAFJAR6cAXd109QrQ2U9Xt1heOiSDCVkHkAI6OgXo6qarV4DOfrq6xfLSIRm8mXUAKaCjU4Cubrp6Bejsp6tbLC8dksGyrANIgWVZB5Aiy7IOICWWZR1AyizLOoAUWZZ1ACmxLM7BhU8Gl156aW/WMSSNjk4Burrp6hWgs5+ubnG9Cp8MgBOzDiAFdHQK0NVNV68Anf10dYvlpUMyKCkpKSnpkDIZlJSUlJRokQxmZx1ACujoFKCrm65eATr76eoWy6tcwrqkpKSkRIuaQUlJSUlJh5TJoKSkpKSkTAYlJSUghKhkHUNJthQiGehYUIUQG4V+1sZPCLGrEGJ81nGkgRBiuhBid/9nba6Zz+jgB93cdPMJk6Rbrnc6E0LsBxwI2MBD2UaTDEKIPYEZwFNCiBellNOllIXvxQ95PQpMzzicRPETwHTgI8CdwJ06XDMAIcRuwLcATwjxuJTy+xq5aff9EZCGW+5GEwkhKlLKmhBiMnAxsAqYA9wgpXw52+gGhp+9hwDHAccC3wP+CfwUuEBKeUeG4Q0Y36sCHAI4wHFSyl9kG1Uy+G7DgG8CVdQ1Wwf4GPBtYHXRvzSFEFsANwLnA/cANwD/kVJ+JfgcZhlfJ+j0/RGQ9ndjrpqJ6grgAmAP4Ezg48B2mQXWAYGTlHINcC+wi5Ty18BbwAvAPCHEkODYDEONRcirD1iESmyP+89NFUJsIYQYFhybYaixCbmtBH4tpfyklPJ24GXgUCnlqqJ+UdZdi17gYSnl/0opXwMuBSwhxLb+l06hrlsdC4DdKfj3R0A3vhtzUzMQQpwKTAH+BFwvpXwu9NwFwOvAlVLKZzMKMTZ1TjdIKRf7f58I/BjVTPcPYJWU8gwhxBD/yzXXhLz+jEoCL/N2rWcM8DDK7SUp5QlF8YKW12yYlHKVEOL3wMwi1ubq3K5EXavbgROklPcJIY4GTgf+LaU8qki1AyHEKcDzUspfBjVx/wassN8fAd36bsxFzUAIcQBwFPAjVJb7hhBih9AhPwPeh8qC4fNye+fSwOnrIafgDnNH4CzgaCHER4vwhVnn9WHgO8B7gd8AdwOHSSmnopLD/kKISUXwgqbXbHv/6dV+p//TwJqMQhwwDdwuAFag2pxNIcRfUHfSBwI7CCG2LkIiEEKMFkLMQjXnXSOEGOrHHa7ZFO77I6Cb3425SAYokcullHej2mMXAF8MnpRS/ht1B/0hIcSnhRBf8f+e58LayOl0ACnlAinlf/2fXwcksH5Gccal3usp4Ewp5SLgO1LKfwJIKZcCvwJGZRPmgGh1zWq+0whgMkDQvFcQGrl9R0p5JXACYEkpDwf+C/wdeDWrQOPgN2+5UspNUDckl/pP9ddqCvr9EdC178auFub6bBX6/UngcAAp5dPAbcBIIcS+ocOvB45HdXiNa/R6WRDTab06J4QQ3wA+CHjpRxudGF7/C4wWQuwrpXwrdPzZKK/53Yk4Oh2Ww+uAHYUQw/NY44nhdiuwoRDiAL8P5O/+cecCI4HXuhRyZFq43eo/ngEc5vd5rBFCDA0dk8vvj4A8fDd2+85mraGsoex1E/CGP1wKYDFqdIMhhKgIIUahes8fBraTUp5Zd36WxHYCEELsJYS4F1XFOzjcDpgTBur1SSHE3Sivg6SUz3cn3FgMqBz6fxuBGnWT16aiuG7vBxBCbCuE+DXwIVQtYVV3wo1FQzcp5et+v9RzwGWoETZIKVf7HeEjUc0sefz+CBgW/iWL78auJAMhxE5CiJ8B3/ELXY//9+DivgzcApzsd1q9gmpeGO5LvQWcLqX8bNChlzUdOI3wn38EOElKeWRenCARr6eAL0gpj8iTF3Tktm7ow/VrKeUVefuy7OQz5j//HOq67Zu3BN7Crae+qU5K+VVgghDiE0KIjYUQH/ObYr+Yp++PAD/OXwAXCiGMLL8bU08GQogPAZeg2vNeQO2+cySozO0fNgL4HSrrzRZCbIaa4LMqOE5K+ULasUalQ6eV/nFPSSn/0+XQW5KQ1zNSylw1eUHHbsHzBCNU8kRCn7HXpJQLuxx6W9q4rZFS9vl3xxuETrsA+AtqtNt6/rG5+f4IEEK8CzWq8HZgCap/6ljI5ruxGzWDnYD5UsrrgSuAN4DPCyHeAyCEOBeV+TZGzfJ8Hvg5ajPn87sQ30DQ0Qn09YLSTWe3m1DNWwgh9gJOA2YCH5RSuplEHY3tgceklFcBFwE3A/sJIXoBhBD/QxevW+LzDIQQVeAtKeXf/N+3R7VpHS+lfFwI8S3UkMR5wPdRF/hsKeUToddYT0r5RqKBdYCOTqCvF5RuDFI3IYQBvCalfCYTgRYIIfZH9a39S0p5m1BreP0V2FNK+YRQQ5dPQ9VmvoPq++jadUssGQghRgPXALuihhR+WUq51K/CfROV4ZeiOoGuAz4KnBWMyBA5nJikoxPo6wWlG4PXrSePTXgA/pf+FcBGqDv7c4BTpJQ3CSHOR7X/n+H3f+yMmlfwFamGMnftuiXZTLQS+CMwDbU8wVQAKeVyKeVZwKnAVVLKfVDLFnywAIVURyfQ1wtKt8HqlstE4LMN8Bcp5aeklLNQTT5f8p+7HugVQnzG91mCahZaAd29bh3VDIQQR6JmZP5LSrlMCDEc6AMOBXYBfiClfKzBeWehOkbOkfka3qWlE+jrBaVb6ZY/fLdgAt8qYCMp5fP+aKEdUCMJT/BrA0cCXwb2B/ZErTt0hJRyWTdjjp0MhBpvvQmqutMHPIGapHK6lPIl/5htUVWdt6SU/xM6dxKqo2QNcGK4LSxLdHQCfb2gdCvdiucWNGUJIaYB+0opRejcs1Bzc3pRa0U90u34YzUT+TI11EYYz0oppwCnoNrynOA4KeX/AXOBzYQQ7xVChMegf0tKOSUvF1JHJ9DXC0o3KN0K5ja77vDdUSOgEEJsAiCl/D6qH2GXLBIBRKwZCDUB4hygBzUmdn3UrNmj/OcrqHa+Q2VoKJcQYgZq3Owo4NMyR+PPdXQCfb2gdKN008JNCHEZcBXwOdTCgHvKHMzxaFsz8Id6zQU2RHXcnItqA5sshNgR+qc+n4NaSCk4byrwddRKltvl6ULq6AT6ekHpVrrp4eb3GRyLqhmsD0zOQyKAaNte9qE6cq4FEEJ8BJiAGu51OTDJ7wS5BfWfMEFKuQA1vX1PKeWf0wm9I3R0An29oHQr3fJHXLd3o75zZwE/lVI+mE3YjYnSZzAXkH5GAzXNeysp5dVAjxDiNKmGPm0BrPEvJFLKP+f4QuroBPp6QelWuuWPOG59UsqnpZRPSCnPyFsigAg1A/nO2W67Af/2fz4GOEEI8RvU6of1HSW5REcn0NcLSrfSLX8MxE3kePe4KM1EQH9bVw01ISJYP/w1YAZqXZAFsmBbyunoBPp6QelWuuWPOG55TQQQIxmg2sfWAV4CthNC/BA1W+40KeW9KcTWDXR0An29oHQrKqVbzok16UwIsRNqYaW/oqaGX5lWYN1CRyfQ1wtKt6JSuuWbODUDgIWo4V4zpZQrUognC3R0An29oHQrKqVbjkl8CeuSkpKSkuLR7T2QS0pKSkpySJkMSkpKSkrKZFBSUlJSUiaDkpKSkhLKZFBSUlJSQpkMSkpKSkook0FJSUlJCfD/Ab13rhzSNoseAAAAAElFTkSuQmCC\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"import pandas as pd\n",
|
|
"\n",
|
|
"paths, vals = zip(*averages)\n",
|
|
"dates = map(lambda p: p.name.split('_')[1], paths)\n",
|
|
"df = pd.DataFrame(\n",
|
|
" {'NDVI': vals}, index=pd.DatetimeIndex(dates)\n",
|
|
")\n",
|
|
"\n",
|
|
"df.plot()"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": null,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": []
|
|
}
|
|
],
|
|
"metadata": {
|
|
"kernelspec": {
|
|
"display_name": "Python 3",
|
|
"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.8.6"
|
|
}
|
|
},
|
|
"nbformat": 4,
|
|
"nbformat_minor": 4
|
|
}
|