mirror of
https://github.com/heyarne/earth-observation-for-journalism.git
synced 2026-05-06 19:13:40 +02:00
670 lines
73 KiB
Text
670 lines
73 KiB
Text
{
|
|
"cells": [
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"# Generalized Spectral Index Calculation\n",
|
|
"\n",
|
|
"Elaborating on the work done in previous sections, this section contains a complete implementation of the calculation of various spectral indicators.\n",
|
|
"\n",
|
|
"It does not contain code to download products from the Open Access Hub[^download_process].\n",
|
|
"It is rather a re-usable notebook that can be re-used for the calculation of indices only.\n",
|
|
"\n",
|
|
"The calculation in this notebook depends on three parameters:\n",
|
|
"\n",
|
|
"- `product_path` is the path to a previously downloaded product\n",
|
|
"- `index_to_calculate` is one of the supported indices; these are\n",
|
|
" - **BSI**: Bare Soil Index\n",
|
|
" - **NBR**: Normalized Burn Ratio\n",
|
|
" - **NDVI**: Normalized Difference Vegetation Index\n",
|
|
" - **NDWI**: normalized Difference Water Index\n",
|
|
"\n",
|
|
"The formulas for the indices below are ported from versions [provided by Sentinel Hub](https://custom-scripts.sentinel-hub.com/custom-scripts/sentinel-2/indexdb/), which implements indices listed in the [Index DataBase](https://www.indexdatabase.de/).\n",
|
|
"\n",
|
|
"When running this notebook interactively, _Kernel → Restart and Run All Cells_ can be used to re-evaluate all cells in this notebook after configuring the pipeline.\n",
|
|
"The path of the output file containing the processed values will be printed below the last cell.\n",
|
|
"\n",
|
|
"[^download_process]: See [](01a-download-process.ipynb) for details"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 1,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"from pathlib import Path\n",
|
|
"\n",
|
|
"product_path = Path('resources/forest_fires/S2A_MSIL2A_20180906T101021_N0208_R022_T33UUT_20180906T131549.zip')\n",
|
|
"index_to_calculate = 'nbr'"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Product Preview\n",
|
|
"\n",
|
|
"The following lines can be uncommented to plot a true color image of the downloaded data:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 2,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"#from sentinel_helpers import scihub_band_paths\n",
|
|
"#import rasterio.plot as rplot\n",
|
|
"#\n",
|
|
"#tci = scihub_band_paths(product_path, 'TCI', '60m')[0]\n",
|
|
"#with r.open(tci) as src:\n",
|
|
"# rplot.show(src)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Define Formulas\n",
|
|
"\n",
|
|
"Formulas are defined as data so that the bands can be substituted with actual values later on. By declaratively expressing the formula calculations computations can be executed lazily and only when needed. \n",
|
|
"This is done so the formulas can be defined independent from actual data, which is needed only much later.\n",
|
|
"\n",
|
|
"### Operators\n",
|
|
"\n",
|
|
"Many spectral indicators can be expressed using band arithmetic with formulas similar to the NDVI:\n",
|
|
"\n",
|
|
"$$\n",
|
|
"\\text{NDVI} = \\frac{\\text{B08} - \\text{B04}}{\\text{B08} + \\text{B04}}\n",
|
|
"$$\n",
|
|
"\n",
|
|
"To define the index calculation formulas in this way, first the basic arithmetic operations `+`, `-`, `*` and `/` are wrapped in functions taking variadic arguments:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 3,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"from functools import reduce\n",
|
|
"\n",
|
|
"def add(nums):\n",
|
|
" return reduce(lambda a, b: a + b, nums)\n",
|
|
"\n",
|
|
"def sub(nums):\n",
|
|
" return reduce(lambda a, b: a - b, nums)\n",
|
|
"\n",
|
|
"def div(nums):\n",
|
|
" return reduce(lambda a, b: a / b, nums)\n",
|
|
"\n",
|
|
"def mul(nums):\n",
|
|
" return reduce(lambda a, b: a * b, nums)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Indices"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"These function are used to define formulas for the selection of indices mentioned in the introduction.\n",
|
|
"These indices are not exhaustive - there are many spectral indices which are not implemented in this notebook. The general shape of theses formulas however allows for enough flexibility to implement other indices as well.\n",
|
|
"\n",
|
|
"The formulas are defined in a lisp-like prefix notation: `(add, 1, 2, 3)` translates to `1 + 2 + 3`.\n",
|
|
"Each element in a formula can be either a function, a string or a tuple. Tuples are delimited using `()`. The first element of these tuples is one of the operations defined above. It is followed by at least one other element, which can be any of the following:\n",
|
|
"\n",
|
|
"- Tuples, allowing the recursive expression of formulas.\n",
|
|
"- Strings, which encode band numbers.\n",
|
|
"- Constants, i.e. integers or floats."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 4,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"indices = {\n",
|
|
" # normalized burn ratio\n",
|
|
" 'nbr': (div, (sub, 'B08', 'B12'), (add, 'B08', 'B12')),\n",
|
|
" # normalized difference in vegetation\n",
|
|
" 'ndvi': (div, (sub, 'B08', 'B04'), (add, 'B08', 'B04')),\n",
|
|
" # normalized difference in water index\n",
|
|
" 'ndwi': (div, (sub, 'B03', 'B08'), (add, 'B03', 'B08')),\n",
|
|
" # bare soil index\n",
|
|
" 'bsi': (div, (sub, (add, 'B11', 'B04'), (add, 'B08', 'B02')),\n",
|
|
" (add, (add, 'B11', 'B04'), (add, 'B08', 'B02')))\n",
|
|
"}"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"An error is thrown if `index_to_calculate` does not mach any of the implemented indices above:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 5,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"supported_indices = ', '.join(indices.keys())\n",
|
|
"assert index_to_calculate in supported_indices, f'Only the following indices are supported: {supported_indices}'"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"### Resolving the Formulas\n",
|
|
"\n",
|
|
"`get_bands` is a function that returns all of the bands referenced in a recursive formula.\n",
|
|
"This is necessary to resolve these references to actual data sets."
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 6,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def get_bands(formula):\n",
|
|
" bands = set()\n",
|
|
" for element in formula:\n",
|
|
" if type(element) == tuple:\n",
|
|
" # recur for sub-formula\n",
|
|
" for band in get_bands(element):\n",
|
|
" bands.add(band)\n",
|
|
" elif type(element) == str:\n",
|
|
" bands.add(element)\n",
|
|
" return bands"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"The resolving process needs a `band_map` in the form of `band_num` → `numpy.array`. By defining the arithmetic operations like above, it can be treated like any other python function - read from the formula and called using `op(args)`:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 7,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def evaluate_formula(band_map, formula):\n",
|
|
" op = formula[0]\n",
|
|
" args = []\n",
|
|
" for element in formula[1:]:\n",
|
|
" if type(element) == tuple:\n",
|
|
" # recur on sub-formula\n",
|
|
" args.append(evaluate_formula(band_map, element))\n",
|
|
" elif type(element) == str:\n",
|
|
" # substitute band number\n",
|
|
" args.append(band_map[element])\n",
|
|
" else:\n",
|
|
" # just append the number\n",
|
|
" args.append(element)\n",
|
|
" return op(args)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"Because the prefix-notation is not commonly used to define mathematical formula, a function is defined that converts prefix a formula from above to infix notation. This should help to avoid errors when transcribing the index formulas, which are usually given in common infix notation:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 8,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"def prefix_to_infix_str(formula):\n",
|
|
" '''\n",
|
|
" Returns a human-readable string of the above data-based notation. Useful\n",
|
|
" for debugging.\n",
|
|
" '''\n",
|
|
" result = ['(']\n",
|
|
" \n",
|
|
" # operation\n",
|
|
" op = formula[0]\n",
|
|
" if op == add:\n",
|
|
" sym = ' + '\n",
|
|
" elif op == sub:\n",
|
|
" sym = ' - '\n",
|
|
" elif op == div:\n",
|
|
" sym = ' / '\n",
|
|
" elif op == mul:\n",
|
|
" sym = ' * '\n",
|
|
" \n",
|
|
" operands = formula[1:]\n",
|
|
" for idx, operand in enumerate(operands, start=1):\n",
|
|
" if type(operand) == tuple:\n",
|
|
" result.append(prefix_to_infix_str(operand))\n",
|
|
" else:\n",
|
|
" result.append(str(operand))\n",
|
|
" if idx < len(operands):\n",
|
|
" result.append(sym)\n",
|
|
" \n",
|
|
" result.append(')')\n",
|
|
" return ''.join(result)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"#### Test Cases"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 9,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"formula = (add, *range(5), (mul, *range(1, 4))) # == (add, 0, 1, 2, 3, 4, (mul, 1, 2, 3))"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 10,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"'(0 + 1 + 2 + 3 + 4 + (1 * 2 * 3))'"
|
|
]
|
|
},
|
|
"execution_count": 10,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"prefix_to_infix_str(formula)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 11,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"16"
|
|
]
|
|
},
|
|
"execution_count": 11,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"evaluate_formula({}, formula)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 12,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"'((B08 - B12) / (B08 + B12))'"
|
|
]
|
|
},
|
|
"execution_count": 12,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"prefix_to_infix_str(indices['nbr'])"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Extraction of Relevant Band File Paths\n",
|
|
"\n",
|
|
"Subsections from here on below contain the actual calculations.\n",
|
|
"They start with the list of bands are referenced by the index formula given by `index_to_calculate`:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 13,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"['B12', 'B08']"
|
|
]
|
|
},
|
|
"execution_count": 13,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"bands = get_bands(indices[index_to_calculate])\n",
|
|
"bands = list(bands)\n",
|
|
"bands"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Resampling\n",
|
|
"\n",
|
|
"Some bands are not available in all resolutions - the band `B08` for example is available only at a resolution of 10m and the band `B12` only up to 20m.\n",
|
|
"\n",
|
|
"The lower-resolution bands are upsampled to the highest resolution band:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 14,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"[PosixPath('zip+file:/home/jovyan/sources/resources/forest_fires/S2A_MSIL2A_20180906T101021_N0208_R022_T33UUT_20180906T131549.zip!/S2A_MSIL2A_20180906T101021_N0208_R022_T33UUT_20180906T131549.SAFE/GRANULE/L2A_T33UUT_A016750_20180906T101022/IMG_DATA/R10m/T33UUT_20180906T101021_B08_10m.jp2'),\n",
|
|
" PosixPath('zip+file:/home/jovyan/sources/resources/forest_fires/S2A_MSIL2A_20180906T101021_N0208_R022_T33UUT_20180906T131549.zip!/S2A_MSIL2A_20180906T101021_N0208_R022_T33UUT_20180906T131549.SAFE/GRANULE/L2A_T33UUT_A016750_20180906T101022/IMG_DATA/R20m/T33UUT_20180906T101021_B12_20m.jp2'),\n",
|
|
" PosixPath('zip+file:/home/jovyan/sources/resources/forest_fires/S2A_MSIL2A_20180906T101021_N0208_R022_T33UUT_20180906T131549.zip!/S2A_MSIL2A_20180906T101021_N0208_R022_T33UUT_20180906T131549.SAFE/GRANULE/L2A_T33UUT_A016750_20180906T101022/IMG_DATA/R60m/T33UUT_20180906T101021_B12_60m.jp2')]"
|
|
]
|
|
},
|
|
"execution_count": 14,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
}
|
|
],
|
|
"source": [
|
|
"# sorting the band path alphabetically sorts the bands first by band number,\n",
|
|
"# then by resolution\n",
|
|
"from sentinel_helpers import scihub_band_paths\n",
|
|
"\n",
|
|
"band_paths = sorted(scihub_band_paths(product_path, bands))\n",
|
|
"band_paths[:3]"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 15,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"highest_resolution_band_paths = []\n",
|
|
"\n",
|
|
"# covered_bands ensures that only the highest resolution band is kept\n",
|
|
"covered_bands = set()\n",
|
|
"for band in band_paths:\n",
|
|
" band_num = band.name.split('_')[-2]\n",
|
|
" if band_num in covered_bands:\n",
|
|
" continue\n",
|
|
" else:\n",
|
|
" highest_resolution_band_paths.append(band)\n",
|
|
" covered_bands.add(band_num)\n",
|
|
" \n",
|
|
"def resolution(band_path):\n",
|
|
" '''\n",
|
|
" Return the resolution encoded in a raster file's path.\n",
|
|
" '''\n",
|
|
" if isinstance(band_path, Path):\n",
|
|
" band_path = band_path.name\n",
|
|
" \n",
|
|
" return int(band_path.split('_')[-1].split('.')[0].replace('m', ''))\n",
|
|
"\n",
|
|
"# using the function above, we parse the highest resolution out of the first band path\n",
|
|
"target_resolution = sorted(resolution(band) for band in highest_resolution_band_paths)[0]"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Raster Cloud Mask Generation\n",
|
|
"\n",
|
|
"Spectral indices can get distorted by highly reflective clouds.\n",
|
|
"To avoid this the products include cloud masks, which contain information about the cloud positions in a product in order to discard them.\n",
|
|
"\n",
|
|
"To construct a cloud mask, it is essential to know the transformation from the coordinates given by the raster's coordinate reference system to the pixel coordinates of the highest resolution raster.\n",
|
|
"\n",
|
|
"These are encoded as metadata in the raster file:"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 16,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"text/plain": [
|
|
"<matplotlib.image.AxesImage at 0x7f5f4cb85ee0>"
|
|
]
|
|
},
|
|
"execution_count": 16,
|
|
"metadata": {},
|
|
"output_type": "execute_result"
|
|
},
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQ4AAAD8CAYAAACGnEoDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAB850lEQVR4nO2dd3wUZf7H37Ozm14JCamQhISEToBAEkFEQEDlUNQRG9bzlFOxnPWKXtWze3ZP79TTOx3bYe/lVHpHekkICQSSEEgvuzu/P2Z3s5vtJZsNv3xer7yyOzvlM88z832+z/f5FkFRFAYwgAEMwBto+prAAAYwgP6HAcExgAEMwGsMCI4BDGAAXmNAcAxgAAPwGgOCYwADGIDXGBAcAxjAALyGtq8JmCFJ0jzgCUAEXpRl+YE+pjSAAQzACYRQ8OOQJEkEdgNzgCpgLXCRLMvb+5TYAAYwAIcIlanKFGCvLMv7ZVnuBN4AFvYxpwEMYABOECpTlQzgoNX3KmBqz50kSboWuBYgqSttUmc5hMeEoQvXcvxwIxmjUtn1/X4mLRzLt/9Yxdm3n84Hf/2KhXfPZvn9X7Lgzll8+NDXnHZVCeuXb6Vgei7V22tISIujq0NPR3MnSUMTqdhYxZhZI1glb2L2dWV88vh3LLxnDsv/8oXl/5yl0/jx9XWMnz+SvasrSRuRTHN9CwAxSdEc3l1L3tShbP5kB6dcMpkvnvnB7hzzb57Bl8+toESawE9f7Sa7KJP6yoaBewqRezrruYv46Lr/BPWeDmw9wugZOXb39OHz6zn7F5MC3k8T7l/Avse/Y9jYVIf39Je3fys4emFDZapyATBXluVrTN8vA6bIsnyjs2PuOf+Pytp3twSLotcoXjSOUOUXytwghPhpRDAa7DZb+AkCWL0/YlwchsbGYDK0g7dtd/jWMjKe34yxpcXuN0EXxucdrzsUHKEyVakCsqy+ZwKHXB1Qvb2mVwn5i1DmF8rcIIT4ORAaADUN6n8xLwchPJxjV5YCoORkWPZpvLik1+k5grdtl/boCodCA4Ax+U6PCxXBsRbIlyQpR5KkMGAx8L6rAxLS4oJCzFeEMr9Q5gZ9yE9wOLjaoXNRFsZpE0CnxTBlFMcL1e3GzTss+3z810dtT63V0nB5qcPzdZ0xGYDmC+xm5xZoYmPRZmYghIejiYpyul8g207ZvNM5n4BdxQ/IsqwHbgA+A3aom+Rtro7p6tAHg5rPCDo/00N/fInjh9MaA23nBILt6yAmxDvcbfjEfShaDUpFFZoft5D5dRcA9deUIujCALhk7/k2x7TPKeKqOx2PhUevawPgtYcfQQgPt/nt4K/LVGqihv1XDUMQBIytrU5vwZO2E7QemjadaFwQOsZRZFn+GPjY0/07mjt7kY3/CDo/01y7NUUgwc2ugeQmaLUo+sC+6L3ddk5tET1eFMPxEw6P3/ZhFtofvsFouu+II610nTaRIR+Vo2RnYtizH+VnTTbHRHy5mfdnjkY5JR3hx01Ad9tl/c7A4RvLiNb8iNLRYXPc369+ij/8eSJCQjxD/7ACo5t787TtxBHDMezeR90vShn8/Eqb3zTR0QgR4Rjqjzk9PiQ0Dl+QNDSxrym4RF/xS394hdt9/OLWQ50PtNCA3m87Q1OT+51cIHfrIRRDt5Axbt6BtrkT/eEaDHv2my5iK4SUrk4MR46irW2i/ueqVlh9yxSE8HCMP+3EEAkl33SvBQjh4WjTUvnzfIn6n5dycFGm5TdNdDQAjRfZ21E8aTtFr8ewex9oRGKrbPvv6NIy4j8Pp/y5DCdHmzi4vUqIomJjVV9TcIlQ5lexsQpNRIRPxwoTRwWYjT16ve38XEm08DMJUTEhgboJsTb7OJtOGHbvI2lzMwDpj6y2aBjpD65gxFXdqyGa8HAMdccw7NrLsXFGMl/8qfvcJmNmzcxu4XTsqlK0mRnetZ3RQPhHa202pS+voO432WjXxjo5yMTP86uEFsbMGtHXFFwilPmNmTUC4zjnFnNXUNa7ND0FBKHcdmDipygWAaR0dNAV7ZlhFYA1W9X/PaZGNtpbxhCULnXaMfLhbg1Hmz3UssuIa7tfen2EwNE5Q8m5Z5Y3t2IHffUhxG82kP6ga801JPw4fMGZURcrXe1dfU3DKXQROvqKnyY62vkSGyo3gxjmcp++hMu26+E70RdwyK+XeWliY9FER9GZl4bmh00OdlB9TsKT4umod2ybcQZtVib6g441lS+Mb4W0H4fXmH1dWV9TcIm+4rf7+WKEjFSX+8y+rixkhQa4absQGOgc8utlXsamJoyNTY6FBli0l9Mvm+DVebUZ6Sz8fIPXfPqtxjFHc0H/JN7L6Jo9Cd2X6/uaxgD6CbS52Ryan07K0yssWpNQPBZNRQ2G2tqTT+NYeM+cvqbgEn3FT/flepcOQuX/GT/Qdn7CX37mVZHegLfc9PsrGPL8GgBOfDQcAKNWQ815eS6PG9A4/p9hzxMl5C9b5fkBIWBT6K8QwsPt/DIAhOKxKGu3uj5YIyImDaKrIMNmeiLowixGU01sLEbT0rLZrqWJiMDY3u4Rv6bFJcS+sQrtsCz0Bw6y56mp5N+wGjEhnpe3fsySvNP5vO21AY0jmAhVfvnLVrnlppnQ+0uurhDQtvPQjdwbeMrPkdAA3AsNoObGqey8L5c//+vvNttPXDDR8jnzSyNoRLpmT2LU/9oQByex4E9nO7/nHtsTv68EQH9ADUzPfUc1+BqbWyh59zbEQc59QgY0jgGosNIsxCEpGI4c7Vs+TiJT/79A/+VQGt7OoPOMRjIWWS2Bm+0QujAYl4+yfhvX7t5PfthRbs9WHcKUsvEIKzZbDhFHjcCwfbdPXr4nnY1jztJpfU3BJUKZnyNu1naRvhYaodx2EHh+PeNTAMJ+Vkfy86tIeyLM9gez70hXp8Wn5v6/XsJN199o4WYtNABqpyap1xBFh9cXkwahzUj3jnN/1TgWJl6mtJ7wbC7XF4iKjyBU+XnFLcgjv6ALIzJKE7JtB/2rbwVdGJroSIzNLShGxWFfGqdNcLrMe9JpHOPnj+xrCi4RyvzGzx+JmOhZPIgmTNfLbOwRym0H/vMTE+IR83MDxAaM04ssn3ty00RGYDh+Qp2iOBkAND9scqj1uEK/FRx7V1f2NQWXCGV+e1dXgsYzo6GnFvpAQenqDOm2A//79sD1o6mdNiRAbEBs6/Zi7cnN04xk1csmeXXNfis40kYk9zUFG4iDk2zyHIQaP2ukjUh2GTLd1wjltgMTvyljfT4+8/4VDPrnSvc7eghlXXcAnK9tl/n0Zvc7WaHfCg5zctZQwa7HhqJM7FYTQ42fNUKZG/QPfo3De8+Jyx/42nbehiD0W8ERakj5JBxh8+6+pjGAICHuP1440XkAccRwr/Z35R0cDPRbwRGTFFoSP/6tDTYOP6HGzxqhzA1Cj5+YNMjme2/wqytL8Wp/Y7tj57KAcXPjONdvBcfh3bV9TcEGZjdgM0KNnzVCmRuEHr+mGba5SwLGT9PtVzF4ufPEwA7hZIUkENwErdZtmEG/FRx5U4e636kPEcr8QpkbhB6/qHdX23wPGD+rl9/Q0GD/e8k4y8eqe8o4cmOZTTiAo2TK1tzcTmecaBXWtjqnhw44gPUO+pOTUKihN/hpoqJcZgd3Bm3OMIxH6yzGQ6F4LJG79/jOz5ugQSvnu9Zzp5J8y36Gx9SxZaIpNCA5GUOtrYZhbjtNRASH3sgh48oai1CyDooDNRfHoWfi4LNBGCIh7RH7rF8nnQPYKZdM7msKLhHK/EKZG/QOP1+EBsCuX6bRPHcMAJpxhRw6NdYvfp1zvTjWSiOJXVFOy69S2Ta/e7m1p9CA7rYToqMQP0uw0WSMo3Ns9tVXH+KCnI388/bHSNrWhTjK85SN/VbjGAhyG0AwcPSGMlKeWkHNzWVkLKzAMNNlgUH/EcA0BtqcYejLD7jeJzMDJT4Gw7ZdDq990mkcoRq2bkYo8wtlbhBa/FKeUtX3xN1d8MsYwA0/jeh5waOeh0ZFcewK56UjXZ3XHEJgzc2d0ADQV1WrQgO8ElgDGkcvwjrpygBOTvRMDC3owjjydi4pC52vkmhTh6CvOWL5fnRpGSnPrEAoHsvexdEMv22VTQi8NnUI+iNHEWNjnbqQiyPzMezYE6C76saAxhFkLLxnTsgKjf7QdoClnGKowab9egy8SlcnaXfZLpX2DCBTFAWlbDwA+x8s5cd7Hld/2LIbQ5x6rCam2x+j4bQcyv9S4jLuxCw0rLl5G7jmDQY0jgH0OnxJoKzNGYbS2BTSMTW+wjzlUPR6jl1ZSn2RkfybVjN0dTSVzYm0P55O3RgtmfevoH3BFCI/30zt28MYvCD4nsknncYx/+YZfU3BJUKZn1/cNKKaD9NJQWZH0H3VI/2+G6/E+TfPQGluDVmh4W/fKnq9ZRoy6J8rGXHbRgC2PjmW6i+GEr2zDl0L7H+glKgDjSgdHaRc4DwiVwgPRykdb8NNLHCdbNhf9FuNY6Agk+/wh5smNhalYBhGnYiw0ruISss53CTUDeW2g+Dx6zpjMuE1zRi3uPYq7TirmOPDdQz52woLN+WUCZbi1v7gpNM4SqQJfXp9d6nW+pqfK/jDzdjUhLDrAOLGXb6fw02Oj1BuOwgeP93n6xwKja4zJiPowtj9ouqz0Z4gEllntOEWCKHhCv1WcPz0Vd9Fogrh4ShNzc530IjsWHfE+e8+wtdC0T3hb9spbW0oY32rPesJ+rJvPUFf86uco6N9znjKz3wRgPjXVxH3bzVa96evdiMmxNu7m5eMswTriXk5LqeL4qgRCJPHuOTg24IzIElSFvAqkAoYgRdkWX5CkqRBwJtANlABSLIsN5iOuRu4GjAAN8my/Jlp+yTgZSAS+BhYJsuyyzlUdlEmDdXe1cgMFJSODgyd6oqJQ5VQMZKVG0vd9sBeN1DZuPxqO0FQ5+cepPj3FGLSIBt7hit+2uyh6Cv6NkNYXz57AMNvV5MAjXp2KVnYuokPXVzGe0OOM3ijYgn9r7qnjM+ue5DLd1+Mdk4Dd33+Hrf98XqnyYRqH4RwbStxy5zbSfzROPTAbbIsjwRKgF9KkjQKuAv4SpblfOAr03dMvy0GRgPzgGckSTKHBz4LXAvkm/7mubt4faWDoKAgwjhtAvrTJ6GrcJARXFE84qcZVxgYMl7WDvGr7XrBJmY8YbvM6IpfT6EhJifb5NwMBgLy7PlR72Xvv4oQkwaR9Uf72JKGLzczeOxRamZ2Lwn/+cpXebJ+GhVb0mk/q5hf/eF6Br3sPJ9IfV0sv89fzl2fvO10H58FhyzLh2VZ3mD63ATsADKAhcArpt1eAc4xfV4IvCHLcocsy+XAXmCKJElpQJwsyytNWsarVsc4RXhM36zxayIiEPNy2HdBGOHVJ1Camh2ul4fHhFF/TSliXJzzk5VXB4ST4CTtvTP0Vds5Q89aH57y0+ZmU740n+dee9LuN6FoNNV3lXm1+uMp+rr93pv+LEJkpP0PgkB4dBh1W1MYeedey+a/nzqNLadEMeKeTUR8uIbEl1e6HAA0dTqmRbRzz+5Fzvfx6w5MkCQpGygCVgNDZFk+DKpwAcwZSjKAg1aHVZm2ZZg+99zu6DrXSpK0TpKkdXFjwhg6Lp38smxGzcwjvTCF4kXjiEuOYeY1pQgagQV3zgJg4d2zAVhw5ywEjcDMa0qJS46heNE40gtTGDUzj/yybIaOS6fo7NEkZsQzfUkxugidZXnL7FgzQz6HxR9+z99KOij4Vx0jHptF4rlTGDe3kNzJWQyfXsC4uYUkzBxD2c8GEa7pstThMJ/jnPvORNCFMffqiegidExfUkxiRjxFZ4/26Z4wGry6p+Jzx9nwMf+fs3QaUfERlC4uIjknyXJPuZOzGDe3kOScJEoXFxEVH2F3T+b/82+e4fc9lVxQ5FE/xRmauLRa4Od33m7XTwulEbz2i8eY+OqswNxTZJjlnvJLs/1/9s4d6/WzZ/7/ifwQYU31lC4uIqUg1XJPCX8+g+HFQ5n67G6mzM223JP+cA0Lbj4FY3u7R/30s88buX3/xZSsd27n8Hs5VpKkGOA74M+yLL8rSdJxWZYTrH5vkGU5UZKkp4GVsiy/Ztr+Eqo9oxK4X5bl2abt04E7ZFle4Oq6l4/6pXJoZ/ALBx3+70jaWsPZc9rLGBQjWzu7WPTBTeTfZMrZYAqFDvvDbG675FEeLJqGobERbUY6+mo1QOr0rS28+9BsBm0+jnGrydAWxNol6YUp9EXbeYpA8Tt2ZSlxBzrRfu2d85kjiPm5GPbsB/znJybEYzgeGBuJdbqAlvOnkv9TeUDazrxk3ivLsZIk6YB3gNdlWX7XtPmIafqB6b/5LqqALKvDM4FDpu2ZDra7RMaoVH+o+4y0c3YwPktVkP5+IotLn7uFpE3dzSjGxXDw12UUHVK478arVTdhjUjLP8PY+7gawPTB708HBY6NT0BMjLcRGoIujM55xb16D33Vdp4iUPwG/XMl2m83BeRcZqEB/vMznHDuOu4trNMFRL+7LiBtJ+bluDXE+yw4JEkSgJeAHbIsP2r10/vA5abPlwPLrbYvliQpXJKkHFQj6BrTdKZJkqQS0zmXWB3jFLu+3+9ul96BINAyr43RTy3l1T8sIOuZrQz6x0pqbi6jc+5kym8ezR8vf43tFY2Ef7wWAG1WOg/mv4UxTs+ep6YS89/1DNrYQOKOZg5ca2sgrbptMmJ772of3radr9GeviKgfdsLmpzf/HrL6dJocMmta7ZaO0UcXeDyNEp0hE1aQ0fwR+M4BbgMOF2SpE2mvzOBB4A5kiTtAeaYviPL8jZABrYDnwK/lGXZ3KvXAy+iGkz3AZ+4u/ikhb7XtfAH+tMnUrlsPJkPrCb2jVWWjEqtU1qZ89D3GEc2c/snF5N7bbcSpT9wkPvmXcTXsx5n/6Ln0cRE01QYj6ahmW03PoOYEI82W035lrq6HV2db0lnPIW3bedtoWJ/EdS+dfCCuMs4Hgh+5v4ONBxyEwQMMydSfr4aLhD1XL3FoC9otXY+H/VFCQhuCnb1W5fzM7SSohiDz333s1N48YyXuOmlX5D5l+7lsGcO/EC2NooL959B04wGBEFBMRgtv9f/vJR1v38WgHuOjOOD16eRs2A/7+d/yk2Hitn0pyIiD7Xx+jvPsXDbZcRf1GCTvUmblYn+oLUN2XcIGoG+aDtP0ev8rBLWCOHhNtnpQV0md+Xm7ZCflzV23fVny3lTiflwkx03d3DETTOuEPYdtIT/n7u9lvVN2Ry6MInDf4vE8E0SqY+toPz+UjK+01NzdTu5t5/AWHOUz1r/dXK5nJ99++lBvZ5h5kS0OcMofK6JRxacR9ZD62x+n/fjDRR8dxX73hgBRgPz/3YB4pDulPcv3vM4X7WJPHc8g1Njd/HJTQ/yfI66Tj4tbjfRFc3sWRrGYDGaL8a8AQYD2rRUmqUSNBNGMeuTwHmTBbvtvEVv8xO0OsTRBYhJg1A6OuymYo6EhiY6msO3lXXz66GpHLx7qoMLOR+13Q0CcV/ssBMa4kj33rqO2k5o7bDJGfLynxfw96wfUY6foKUtnNTHViCEh/Pw+a/QFa1BuyaWw3+LRBiWaXcuyzn7q8YR7LD6P5av5dI3byLnLtel+8wJWG7au5OniootUxmheCxifRMHHo4mIaqNQZGt7P0il83XP8nurk4Wb7yajMX7+LRcXZ0Z88RS7rhS5sLYw/yxdiKfVxfS8NNgcu+0vb52WBb6AwfteAzAPRwl+wX75DxiXg4Zrx/lgqQ1PH7KTAxHjmL8Kgvhd0mIW/ZhbGpyWPHdkTYjhIfTOWMsus9tBx5P4CgxlKALg3H5KOu3eXWuQ7eXkf7Qiu5zakT0Myeg/Wo9LedNJW5rHdQe49P6F04ujcO8Ph4sfNI4nvBj7r39zPaAA+/+2iajtLJ2K/r9FWReuIeYhdV0zqpDo4dpd93Auf++lU1TXsPwsaqhGBQjbWPaWBJXR7igw6gI1B6OZ+7pGywjXds5U9R9UxLUC1iPgBoRMdl5DdFgt523CBY/R0ID4NA14xEHJ1m+T3x7L49lfMUZUV0M/m87C++ezRcjP2D8k5upuUy1KfQUGvVXlzqcZrTPHkdDnm8OZI4SQ1XdNhm2dGf+8rTt0h9aYXtOowHtV+qydfQ7qzHs3ue4ZIMJAxqHBzBOL2LCE5ssaekDDTEujqMXjmbd758l9+1f8MyZL7Ns/YXsmv4qANs62xgdFskpWxYRM89kNS8ZB6u2AGreBkOUEU2HhuG3r0SbPZRDZ2WS8vQKtKlDMLa0qoFper267j8uz3JszziR/g53IfueQDtM9Rowa3LPHfiBu6t+xs9Tv+PRaXPQH67hop2H+P2355Dwk9aSl9QaYmKi8xfPjT1EHJyEoa7eI66W+zXZbZxpUb7ipAurN3vmBQPaE21sm5Po0b7iyHzEhHiv+Cl6PYNfXs9thyeS/j+49/dXEf+xmjquSzFQY4hhS2c7Oo1RVU01ouXFB9iz5FmeO/MfxFaofayvqCTl6RWIg5M4Pj0bY3OzRRNSDAZm3jqK1kXqnFyIivKr8npvwJ++DUQgoL6yymb6d82Vyzg+o5FHTz9T9cK8cxbytLGMuG6NRWiYkwWb4Wq0diQ0rNMkOnQn7wGzcDPfryY8nMp7y1jy37GWKFjzPr2Bfqtx9NWqSv3VpSS95NzOcfjWMjI/qcW4a6/XKwPmSFtxcBKGY8cp//MUBm1TGLS2DiUyDHZXIGRnIjS2oK/qjnOpvLeM3JcPWh5289y67oMRJN0fSXtKOFEfbWLX3yagje9k+KWbMZSNR7txDx2lhZwYHkbyukav58mOEAgNptdXVbxcAemJnvzEwUkIkZEYUuJ9akPNhFFU3Qvp53puAL92935eGJFrs02bPRT+aUA/U3029j9QSu5dK6n4cynZv3Zhm3NRkuGk0zhOu8p5GvnegCY2FnFICvEXVdut87ecPxVtphpek/boCoSGRrf8NNHRdiO9OTzfUFcPRgM5d68k/rVVGHbtxbhpO8bWVsa+vpszP99ic9zQ36+wGSGPvav6CKT8WsO9r/2T5U8+TvsZ40GBPae9zMylp6H5YRNTVhzn2qffJeWVjYh1jU5XARouL7UP5HO2YqAR0aYOcXnv7tDbfauJ9C+viR2/IYMxHmtwKjQarih1mXhZqDpC1u1tXnF49DcX223TV1Ry6pfdkcK5d61ECA/ntkXLbfpL0GoRRwxn/wOlaKKiqPxtKUduKvPK0a/fCo71ywOXD8ItBAHj6ByUtMF0PpOGPlkVItrsoQhaLc8//DhXff09HWepruJduaku+XXOncyuZwr465svek1l88/H8MLfe4Tx9HiJ279SDaPGLTv5pHE8Nx48k/CP1jK/eAvbOtvY+H0Vgi6M3ydvY2rEQfa/MgJj3THE+Di0aak2Fb3EwUm0pAsM/iYSBMES7Vt/dYlFgO5+vtgiWHY8MAylw96Ip80eyp6np6KUjVeFpgusX77Vfd1TB/ftKXyt6mZGz741bN9tswrTE7rFrpM6GerqMewt94pDjOw4LP7rOvsYmOVzimw0Ck12Fhw7zv8ufghjWxtPXf48jy97DjF5sMfX77eCo2B6rvudAgVF4cD8GDoGRxL17mqe//dTKE3NXPzZjwz6LpYF397AH5+6lMiDpqXXFZtd8nvi+acYmnaMZTfd6D2VdT+R+pitMU7Q6my+pz3a/fvGs7J4adgXABxYEM+vzr6SvOGR7Hp6PEcNLcx6/zZ0G2NAo0EZlk782x1kvXwQMSEeQRdG2dfVSIu/ZcWGAuqvKqFrvCosdC0K9332BggCm876G0pnJ2JeDlnva9jx1zy0GelU3lfW/XJrNCSt11B1erTLlwwg/Z5ZHLhlggeN4eN0xt1xbtyt7frWzflizjwQtFIZxbW2gkPp6LCZ1gIY9pZjqKvnnN/cDorCzz++hjt//wv0h2s8vk6/FRzV2z2/Sb8gCGhzs9FHGy3Zutd2ZND4birnx9Tw4rDPUBR4ZNnztDzUvfzmjN/RG8q4/KFbOdoYQ+QnGxC0Wo7cVEbzBd0ORPv/PQFQPf7M2atdwfqhFHRhNiOx/tBhpj6wTP1ccwTjTzup3l7DqAdqueTiGxi8TkPGX1fSOWUELdmxtBu0fL5hLEp7B0pXJ5lhx9jRnMo/5v+dE/kQVlFLs1RCzalGpoTraL5gKkZFAUXhvA9WYtQKRBzSUXf6ML646kH0MyeiiYoi5d/1xFR3IbZD5e/KXN7PiZV1nLFojVvX716DG/uH189eECOfveGW8Kpq98i/cbXls6fot4IjIc1FgpwAQpuRDkYjH1/wCBqTtfvujy/ix3Hvcsp9NyEKAuVzX2JWpIFlOV+x75ESDt9a5pRf+nsVpDyzkqzzf0IIC+PQsil8+KsH+cMDpmmLILB9xksA7Lo2nvFPb7FX2wUB44wihwmTFYOBfQ9azcEVhSFP2mooCWlx6MsPIK7ZzqDX1UA83YkOoqpa6Lg8iuFv6C3W+j98togDT43gyi+vJv+VOoy1dXz86GPcO/O/qr/JpQ2c+uSvAHjt5rOJOdBC7stVnFjQTKY2hhM5YRy9dDxdRpGmoTraJraS/d9jHF9S6rTNE2uPsPWO8Rh273O6T18iWM+eM2gz0p1qRcHi1m9XVZZOu13Zs6Ki16+T9GMixzsj6fptCuLq7ZbRPXlFAlV/yufbF/9us//UO6+nfZDA9LWH8ISfGBeHfkwOYksnxs07AKi+q4yMB1YgaLUcvmEKjWM7GXF1t6ehmJ/L/j9H03Uomryb7ee6ruIgtMOyyMkQWTkmjegjBiKXr7H4HFiXHXSFQ78qY/XNj1P08jISdyrEv+Zgvi0IHHhzDBnPhRG2agdCRiqE6fj4izeZ/LvriTxmpDlVJP29/XYqcn5Ztsu2s85B0Rdwx683oYmKomrpBLI+qnVY8jHQ3E66VZWO5uDMGX+Z9hVL0lei+WETSlenZY28bkYr4R+vpUtR1dAOpYvPW3VE1umJP6Cno92IUGxaNbGaOhhOm2hzfkNjI8KKzRahse+REjIeMHn16fWkPbuewpt22BzTlR7PoHejybt1rR3fIzeWoT9Y5dS4WD8tg9YR2Qx+cQ0xP6ojuqGhgSM3lVmEhiY21mWbpD+8grk33kj2b1YS/297Dip5hWEX/oT26/UYW1sx7NlP3iv7+UtdAR/f+zBhx/XEVuvZ/Ugqh24vs/F0ddW3mthYNCmeG/F6A0F59gTB1rfDtOLRPG8sZ1y0CuVAt93Cuq+D9V70W8GRNNQzhyyfYXrZ71l2Ha+e2r1qUPuK+tAKogZxSAplv7mBqyuncUPVady29QIif9hJ7PpD6O8agyIKCFotB94cw54npyLowqid4HopMKratkuUjg47Y6Lmu43EvrHKMnc2L392zZ5EbJVpPt0jD6l5JSP+9VWMvrYBcWSejXdiZG13JO++X49RV1BMQtIhz/dMGc+s5u925Rt6aLN7rs7j36/PYupHt1AzNZzI5Wt4qeRlllz+GUcW5VmOt+5bc0JnITycyt+Vkf6FgqG6BkGrpfKtsXaOV57Cn1ykvf7sAShKt+1KI1qEeuzXO9nw60lo4mLVFaqICJTObmERFG7046mKlHWN0qsp6q2chAStluqbp5D+8Aq1pmlrG7t/NRxBgbDhjcS/FUvChlqUqsMYW1s5fGsZ4348Tu1q1aHns0Ob6FC6uGT/fLbVpCGuj2XYm9UY647ZxLP4isQfB9F0aSzGuCiL5mIN/emT4J5atLPVDOHJUwqpXeOiOpgLhyBnEAcnses3+cTt1pDyjL0Lts3prYK/xBHDqVw0xKJlaSIiiE8Kp+FwM2JyEjXnDif5uZVQMo5950VjiNMz4hdr2f3PSdAmUnjrFos9xjo9o1OeQ1Iw1Naz52+Tyb9htVf3aEZiRnyflkeoemc0mRfspOXcyWjbjIR/vJbjS0pJ+NcqBo/LoX57lUXoeOO+7ggn3VRlzKwR7nfyB1YjqaLXk/WKOp80HqlFiIrEEG1g2Mcd3D76C065YzXHHhfQxKuGqYT9egoLun0Vjhpa+KlT4e3hX/JYkcxPNz3Dwo/XUvQ/zx4+V445gi6MN3K+5s6vljsUGrXXlfLiP5/gtRH/thjUCgvtpyKCLoyOM00pC01Co/VcB6HiVmhdNNWSkGbXb/Ip+N12IhqMLo8BbIK/DLv3WYQGgLGjg2H3zaLzs0zO+fYn0v5ris1ZtYURf9lJwQ1q2cnhryiMvGeXjYt50z/CnU+zTBpkV2EGgkag8Le+F1Xq9WfPDTIv2AlGA9HvrLZkmRuzdCt7npzCor9HUn7vJMQhKYhDUugaPZR9j5Sw+7kp3l3EjY9MvxUcq+RNvXp+61wa0B1JaWxtRV9+gIKbNiB+u4G3Zk1m26wEGn9MofLSXCr+WErYcb2FX8PlpSzZI3HpKzfTbGxnVqRq1Ls2/hDv/Xeaw2trc4YB3fELzedMcpoxSunqJOf9a7n8u6sd/h5bbeCBI3N49UQRWpNtwFHbKQYDkV/bOjYpLtwZmhaXIBgUDMmqyl/wwjEMjY2WIkBAt43HQzRdWIJSOo5Rp7zAN6OXMymigqPzc1FOmQCothjLSPrtBpuEv3uenkppSjlKuyqUtJkZ1F9TijYrk/qrSxHC1LbU/XQARa/3y7ja28+eWzhY3q05JxpFVDgR/T3Zv1uD4ahpgNNpWH/hY3x35qMOTmQPc1TwsSvdeD57zzo0MPs6174A/qJhVi7aLOeJTMxzTn31IQwNDeiKG9h6yzMIBc20Dglj9u/OBiB+SRUaqZ2h962g+IVbKb3vBvK+uZLRKy8h57Gf0Iwp5MQltp00TD5C+9lT2PPgRDrOLKb6dAUl2uS5mRCPcYZtAaIR169lxJWOM3lHfLCGqgsG8925Y9DXqB6MDtvOaLALEIt+27EqrxlTyF/+/AJn/ekblHU/Aar3pF0b9az25saxKv6DLWjWbOP4G1eQ+84v+M3CJRw7vR3FTRo7gEGbNNR1xFgEiz4zicSLqtBXVZPy310WN3NzHI27zFqu3Oa9ffbMAYW9mYRaX3OEkb/ZS8Pnv1AFi6KgLz8AwIJfLuO092+z7Cvowiw2u57JgYwn1KnzQ79+3rWm219tHL0eVu/NPF8QeOfgSk5ZdyXN5fG8uvAZrpR/Se5bjYjVdZYX1hG0aakcWJJL1qcNGDfvQJuRzt7rhpHzpw3senosKf/TkfDaGlCMHHxrNNcUruCiuC1cMdSxthIMGKcX0Z4cRmRNO8IK9xXrNdHR6IvyqZwXSfZvvHM0EpMGYWg4YTPKdswvJvwT+9WcniHl5uVl47QJhO2rceoZqc0eitLWjuFId1kBzZhCjD+5rhLvKfSnT0L79XqEotEoG23jWYSi0Rw+Ld7OG9gOngbmuXluGy4vpSVDsKS9rPtFKYOft+8TbW42+v0VJ5+Nw1xYptdglZPSHcT8XCa+egvJT0Ry69yPmBxuYOngbJ5993n2PD7EdqQRBJvlM/3hGsKPKzzx/otU3lvGz7/9H4N2KAiFuWij9CRua7SMICn/iOSNA5M44293WI73pRC1v22n+X4jUe+u9khoANReNI4DZ0cy7JPuQC5XFe6s+Rnqj9m9MFEHHNuGeuahMGuFNWVRrgcBg9EylTHDldDwtv3MdV16Cg2AuolxzL7UeTnGbkI9hIYTG4S7RD76KBj6+CbLd0dCA0C/v8LleQY0Dhc4flkp427YQuVU17EV1tCMKeSTz99gVbuBTG0brx6fxCvvn07u2ycwbtpO7fsFxD8TR9ina2ldNJWod1dT/e5ofip53XKOVmMnU566magaRS3XZ40eI4omNtZ+ZcaHVRHXN+VfGLo2Z5ha89XESTNhFA2j4xj0wXa17oyXEEfmO3R+6o+ov7qUrjjBvcbRA8FKwDSgcfiA+gkKk2IrQCN6HHJs/GknI19Yysa3/sySq5bx/WUT2XXVszy1/AUATjRGcbRIx/HLSolbfZCaW8rIeEBDpb7Zco4oTRjhp9RxfARU/LGUo0ut5tRWAkEoGo0mPo69j5eohkCTAdVd6n1v2k7QakFxv1LiCvryA7a8O7rQRwhODZTu+Bl27nX5e28jkM9eyjvbSX9+k8PfXCUndiY0gvFewIDG4RLa7KEoTc0oaSkYt+3yehRvuKKUh3/3LKdGwP/a4c+5EwC1BOQ3V0yhaXgsXHWUTr2W2L/F8cU/X0AUNKzp6OLGe28iaXUtr36p1u++JOsUu/Pvf7CU+xbKHOgYTFn0HpY9eR2pj3s3crmDZvxIh8u8/sA4bQJHi6Mcj7IB1JYcJQt2BYfamw8I1HmMM4oQf9zqc10bTWwsmrhYG9+WPU+UMPKBCsf2Hgdtf9JpHOYCwb0JfUUlhvpj6nzXy4d5ztJpJL21hWveuJ6zd8/npseWWn5rN+roeKCZhO8rEJ8bTMyTcYT/sJ0xz9/AxHUXct8ZF5Lw6koMu/Zy6nO3c+pzt1uONS9NAsRUCiSJzfzv8snEatp565aH1PV7Uz4N86qQJjrakivEzM1TeCw0vMiNoflhE6mPO55bi/m5Nvxc2ULcQUxMoH2Brf+CoAtzmg9EsLoHVytq7tqv6S3niaK9QdiuQ14LDWtuxqYm9IcO2/yuxOqpuMJJygcvnvF+q3EsTLxMaT3hf37JgKGHtI6Kj8DMT9CFoei7uuf4sbHsfGgkGdl1fDTmdaTMUps09d7YE5Sy8egq67jwy9X847ZziVqxl7Ypw4msOE5bdgJhn62j9dypVM/C4ilpzS1QOHJjmU0Urj+jrpmfOLqA42MSiX3TA+OhE1jzaF00lazbdnP8qiQMu1xPdxyVIujJzxmUsvEODceuzukrxLwcmyRA7rgd/HUZsQcVj8LoxcFJfHr0uZNL4xg/f2SfXl8/a5Lls1iQh3HaBJvfrfkpXZ02QsXY1MQp43aTHXeM0/6irq9rok0JaqeMtuyniYqyjLjOPCKFFZtROjp57t7zCf9kA4aGBiK++4mj05IJ+0yNqI16bzVKuNGirdi1nY+ZtKwx5Glbnw+zI5YvsPCrriFxVbXrnc1w4iNiLbw++dsTrNowAuVQ9/K4tQZnDVcvuIWfk2s6W22y9HEAYdhfafPd3Xsx9OH1JDoLTOwBJd255tRvBcfe1ZXud+pFhK20UuGNRru6Gu741c1so25GKynPqS+c2Qtyr6Qu1R67spSGc8dhMD349eeNQSzIsz2J6cE11Naqo7JJUzG2t9slVC5YupnWtHAbboJWy7ErS6m/WnVAMxcl9gk9tCR/RlYzP8PxE2ouVTeOYwCNF9o6V/V0ttKmDmHcOzeTf+NqG2Hyyn+e8ihZkiN+7uqr9oS1p2vA0KPd3T13SkeHOv3xYLBwVQaz3wqOtBGBmUf6CusVAcOe/XbzQ3f8lA41w5ZNdGl0NJHZTWgzM1j1p6eZcst6Os6cDEBDIUx8cxcANTerqyzioASn5y9/wJQoRyOizUhH6eq0eIKmjUhGm5aKEBbGubd8zf/uewKAthSds9MFFT3bTtC5XtESwsOJ39vSHZqvEYn43jZjuL7mCPnL7Kc8y5sLmPWCdwZlM79gF+MG3L7wnrwXwuQx7Hlqitv0Ca7geVrjEENzvee+FX0BX/gZW1p4ZsLr/PBxAaKg4bNPJjN8z1GU6GiG/2Ytv79sM2fpStA1q0LKWdSjJiqKRJOvkaDTorSqjlfajHS6hiXTXH8M/eEaxORkytsE/tmopugb9F0l5leh/ppSEne20ZoeQXR1uyUDezDQs+3crYyIiQno125FMTvWGQ1u85qa8cHcIlpHpRKGbUlGVysyffrsCRpQnNvAPOG2T4ql8A/7IMz3gaLfahy9DgeSvVkqQZs91GmeCmHSaLRpqTbbLBGnHuLqt6/nP6/OYuTzS8n+7UoMu/ex7zfjEMaMYN7OhYhZ6SS96NywJSbEgyCQ8K+V7P7HZJSODiZ+U8fRX5ZRO2cYYnMHmng1MM1QW8ve343kowtUDca8bKfNHkrKd0fYd2E4TVmaXhMahpkT3e/kAcwu/a4C16xLOtoce7DKYguyRsZ3OrtAx5BAAPKX5t6xEsORo24dyFz5LvmtcUiSJALrgGpZls+WJGkQ8CaQDVQAkizLDaZ97wauBgzATbIsf2baPgl4GYgEPgaWybLscrknJsl1in2/4WC1KWHNIboyBqFok9B8Z9/odRPjGGyq4xmTFI2gC+PQqVpyP1djJsT8XHVa4wL5L9fbBYzl3rcBwsLQzGlFbzQgxsVhHJ6FsnEbex8rYcRvtmJsaaHjrGIG3VlBy6nqXDpuSxh7Hy+hen8TXakK2R+1sOeuCM4c3s5eU8R82Gfr6Pko6g8cBEVB05FCU76edAcxFoFA2PYqu2tD7/StEBsDVhqaGBfn0mt1X+NgmhekkfTiUbvfevXZ89NLtye3/X8tVQuV+3BeVxXlAqFxLAOsF/vvAr6SZTkf+Mr0HUmSRgGLgdHAPOAZk9ABeBa4Fsg3/c1zd9HDuwNXH9NT6CsqEX7chG57lUONJOnFVRaj4OHdtTAun5y7VqLJVcPkCdNZCjfZSHMr459h+2679IJKR4dq0DMaEIrH0jV+OAfnxXPgD6V8c97D7Hx0FJX3ljHqvq389GO3ATX18RWM+OcJEl6PIef3a2HVFvKu3EGb9pDdKpDtBVWhOWTMUSIOazk2rtuXovw/49Xs4x4YLN3BOqjMGr3Rt/ryAzRd2B2FvOsPo1zuH3HWIaeaXa8+e468dL1Y9erJLW2FKiyUUu/LfLpaTvdLcEiSlAmcBVhXFloIvGL6/ApwjtX2N2RZ7pBluRzYC0yRJCkNiJNleaVJy3jV6hinyJvq2q26NyAUjabx4hJapubY1TIBbLSUvKlDLZW9zNm6Ddt2caIkk7ZzptB16niMM4qo+FMpDUtsnZTOfOpbBK2WlvNtE+looqI4XhDNH195kQ03PEHsxHpOl39F4bMtDPvwBJ+sHYdxaDvarEzK/6IaR42bdxD9zuru2rGdnRzYfxrG8O4XX5g8xuH9xswvJ/vRrQxep2pX2swMHit+k/ahCRy/1MvEMF6gt/o24ctuTS7vFtfZv1ytCrniJ+b7We/HkV+VF75WPblFLl8DYDfdFHRhDrPkewp/NY7HgTsAazE5RJblwwCm/+aJYgZw0Gq/KtO2DNPnnttdYvMngXWD9ggaODZKIOLDNW6XG53xi357NVGfbEYfJRK27yjz568l+btul2AxL4cVDbmIWRnEfma7MlD0YzPv3/8IXzePQovIoJ/tY9iEQ+jjw1HWbyP9G4Hhl21h2TefYRja3q0VCIJq7CsdD4rCoWs+t1lOVjQCQtFo7KAoGJuaMGxTV3OMDcd56vxFtA7R0ZgtWGwAHlVd8wK91bc2c3o/HB9d8VNqgq8JW8Oam7PIbnFkPkfezmXKxxU+X8dnwSFJ0tnAUVmWHWeQsYcjfUtxsd3RNa+VJGmdJEnrJl5bwNBx6eSXZTNqZh7phSkULxpHXHIMM68pRdAIlqrn5lDjBXfOQtAIzLymlLjkGIoXjSO9MIVRM/PIL8tm6Lh0is4eTWJGPNOXFKOL0DH/5hnqOe6Zg7J+G8taYwDVtTcmYzCli4tIK85n3NxCcidnkTs5i3FzC5l9/TRKFxcRFR9hcQM2ByD97LZTifhwDXPOz2P/GdGUTk8hMSOeorNHkxHVQdo98YzI1pGaEWlzT2/vmMhrnz7Jj3NzePyrh8BoYF71KYjfb2LmNaWkf7uN4nNGU6AkMf1oOgUXnULOaaOZ+IvTGZSVyPQzhqGL0HHJI+dgbG1V+QgC55w3grgna5izdBpR8RGULi4iOSfJ7p6SUiKYWhhOxwI9iyamoqQNZuE9czC2tVnubf7NM9BF6Ji+pNhyT97202WPL0IQNZ73U2SYTT9Z/3d3T8k5SU77ydk9nfObuU7vacaFYxze08/uOcOvZ6/nPXVcN93hPS28Z47lniIjBOb8crrdOYx7D3DF5gLWXzPRcT/NGUlmST7Fi8Y5eg0BP1zOJUm6H7gM0AMRQBzwLlAMnCbL8mHTNORbWZYLTIZRZFm+33T8Z8B9qAbUb2RZLjRtv8h0/C9cXT8YQW7uIEwajaa5A2PFQa+CqSzHa7UIYWF2qwHatFQ1CMnbgC/T/jU3l5H6+ArEEcO5+sMv+LRhLAdP6VCX8iYUgN6IsnEb2rRUjp2WzTn3fEWXIvL9OM9ye1TeW0bO2/UWTcRRgppgwttgtpMBYkK8a4cy07NgnYJAExFhl+XN4bkHJ3FiZh4xb60OfJCbLMt3y7KcKctyNqrR82tZli8F3gcuN+12ObDc9Pl9YLEkSeGSJOWgGkHXmKYzTZIklUiSJABLrI5ximCFD7uCsn4bhj3lDh9aT/idkCaz60Fbo5UmIoJjp2UjJibStrDYEpDlibOOYCqJYI6Q3XnjYOZH1XFawk7VfT19CMr67Zz7aC6V95ZhTEqgMVvDK/IcfjzFgeOQEwPo0N+vsAgNAEOc80rsvsDbvg220AiFZ8+Z0LBwMw041nlLeiYrsob+9G6vYUNdPTFvubYB9YYfxwPAHEmS9gBzTN+RZXkbIAPbgU+BX8qybF4fuh7VwLoX2Ad84u4iy//yReCZ+wInS1ye8Ivf0UTilu4uEOPi0KQNoTFbg5AYT1O6Fgzq+RvfSnZpSzBOL6LuClufkYJ7tjP6gxtIEpsxNDaiP3AQzah8Xm7fDoKaOyRxRg3Z8lEMY3NBIzqNHHVVY0Xz3Ua39+oNgtm37Wd7b+QNmWfPAVxx67n8rBlTaHHNF3potq6ig03797nG7xOe+s0LSih34MJ75nj9gAnh4RgnFSIYFYT1O9n94hjyr9iIIIqckCYT92/nUaL7HiplxB+32T0cYkEeQmu7pSRk57xiLpiUqHJTFI7+soyUZ1by630buWf3Ihq+SyXzfnsX7EBGdooJ8SjZGQjVtXbp/sC3tvMFyikT2HdeBHm3ehd9Gyx+vsAbbpW/K0PsxKY8hRniqBEYtu92OlXpt4IjFGwcnkDQhaEYDF4535hfUmHyGHZdF8HIX1c69XkwQzssS3Xc8hEH5LG8N+V5rrvxZiI+WONwHzEuDmNrq98xGpX3lvH2lY9w0aaryFxSjaGxMWDJb7xBX1yzv+GkS+RjtjiHKsz8jl08yZKa3w5OHHvMI7uy7idG51WjNNvHH2iiomyW2zwRGmanM0dtl3Plfm495xoboSGE2xY4MjQ2WopOgTp9EXRhduUa3CHniW0suf9WMq84hBAbw6FflbHz8QJLBftg9a2vQiPUnj1rd/pgceu3GseZURcrXe1dfU3DKXQROrrau5yujIijC1TXcjftr80ZZqmPIUwabXEq2//vCURsjCLr+Z/UFzoiAsVgdDqdOHFpCUdOMTLi+jXd3FxATEzkxOwRHD67i4Klu+zr10ZFcfDGCeiaoWMQZP1ppVflJJounErsG6uouaWMDb96impDK9ftk6j4Kpvch9e65deX8KT9gglrbTPQ3E46jaNEmtDXFFzCws/Jy1Q7ZRBdsya6TY1nFhoAwo5ymi4soWN+Me+UPUdzrp6uIjWyteJf+YjpzosItaRqEJtFW25gk1pPzMvpPiAlidjrq/hp1nMcvmq87QqLINDy3hCGvXGQyHojWX9aaXGl9wTCpNFE1HVR94tSBp1ZjShoyBCjMJ5Rx7UXfdx/+jZEYK1tBotbvw2r/+kr32t/BgPu+CVtaqT+j12En++5wdHY2opRCxE1rSz89CZ0CR3oDjdz7NISBsUeRX/gINrUIey+NZfcO2zjLNIf6dYIzNy0mRlE3lKNuMFUrPl4t2FVPziGrgcTmXDFNQz/7hhG61q6peMwvBiB/sBqYg8eomvWRI6nhZHwryo8gbLuJ3RA6o50KmOzydn/c2J360jrWsFnF5WQJDaw/98TGL7kJ4s9RRMV5VfZxkAilJ+9YHHrt4IjuyjT/4rhga4/YgV3/IS9lYT/cxTGdu86Ov71VShA4a3qi2QAGqUhrBvzBou0ZXTlpHLp3O9Y+dtYW/8GRbHcr5mb0txM9YkM0obpEKoP2eT3EH7cRBiQ82l3PIHFaLtiMzGoWcpOjIC8kgOkhLXT9L8slKgI1XfAg7bVVx8i/eFDNpGqxi07yVwwhg3bo2haNJno6naqT42iNbuL5BVakt7ajCY+zmlVtmAgIM9eLyFY3PrtVKW+ssHvcwQ6xsIa7vgZm5qIfse1k43L460zkEUqTF17BYpeT8ojB/ioajQnzus2WGqiotRMYEMz0WakW7gZjp8g/bxdHldk62k/EYywbOGHzE7eyZ8z3+fWrz+mM1U1prYt9DwPSc8l5PoDx8j60wpi39uAuHE3rdldDPle5P7fvUD5y8NJfMex96MvVe18QSCevd5CsLj1W8ERHuO/t6KnWaJ8QSD4QY/IVSeenNm/Xkn6eaonZ4MUzeCLjtj4fNS/lYEYH8f+y7Mo+uggkUlxIAjsf7AUYaLz5Lbi6AKnv7WcN5X0q/bz0Sl5fDkti+G6GKZFtKP9fgsIAjErK7y7USuY207RdyGEhzPiF2uJf20VP//+CnRrYlm5upDGi0vsVqWMQfIgDVTf9gaCxa3fCg5deGjPsqz5eTMSNl9gG0qPQemuk5Lm3Php9hPRV1Xbj+C7k+gaPQx9jMLHz01DFx1B+V9K+Nms1Yh1zpPZ7JccZ80CiH53DXvrB2NoaMBw/ATn7JlL2R9uUm0SimLxO9FmZlhqmWiHZbm6dQssbacoGBq6R9D8K9aT/vAKTi/byt2/f5Xqd3rk1AjSCmEoP3vB4ha6LeAGxw97X3M0mLDmJ0RHgQfBRQBxn263yVGgbNyGvmw8AthU5PIGebetBaOB3B8FTlw8lYbqBnSNUdyZ8j2XVZ7q8BhBq2XY7+2nUpZAKUUhU9plCWNun1XPYP2RHjuL7PprChFbhjH0w3rKFyWRsCeD+I+2OfShEEcMx7B7n9u+/XJ7IV9pCtCGe+aIZhcQ5qdtK5SfvWBx67caR8aoVPc79SGs+Rnqj6mjrgeah6MXSrffT0OgeUVEUeCSOuIuG4+2DUreu9XOo7XtnCl0zitm1/MTHHq7alK783Bae5A68ibVhOmI/18EeWfuo+KcJP5x5ZN0xAsobW1oMzPQZmZ0ZyYHS1pFd3074sr1xGyJYPjP1UJEruJoAFpPGWHzXdDqME4vQjPet9o8ofzsBYtbv3UAO2/IlUpjbbP7HR3BbCsIQOJXZ4hLiaXxaLcQ0MTGQleXR2HNvQlBF0ZsYng3N0FAM7bAUkOj/uelJK9soG1YLOEfeVa4xxX2PVxCdN4JMm9p5dgzWuLOqkApHYt2634uWbuN+96TGPHkAQypiRbntri0eBoP997KQPT/kvlb9nuc9cgdpD7hIC5Hq3XpVh+XHIPPz14vI9DcTjoHsEkLvc+haIHR4HUxHW8x6We26fiMTU1QkONkb+8h6HwzgildnTbc9rxcRNU8dcQ2nDaRpL+vRNlTbi80rA2RXuQbzbtrHRlX1qCvqCTuzP1q2/+4CUNTE/dvm8eeJc9y4MlEyu/oPueks1znA/UXTV0RnPbWr0h9YgWaCaMsNiQzFKPrwdSvZy9AcJbuMVjc+q3GcYZWUtx1cF9C0AhuH8DAXEhQIxmt8mNYQxycZFd/xZqbWJBnqaOqmTAK46budIVmV2ZtzjAMB6sto7D1Mf6g4k+lfHDZw3zYNJb3D40j/IwKlZ+bEd9fWJ+/48xiwk50eVUCImh96wI9a8aaEWhuJ53Gcfbtp/c1BZcICj+NCIqCcNxJsJZGdFi0yZqbtQCwFhoANfMyOXpDGcYjtYhW9WIMu/aiiY6mbaF/CYvzni7nwr/ezge/mkVdc7SlVu5Zn53ps0blCayFUvjHa50KDWc5O0Ph2XMkNCB43PqtxtFfwuo9gd/h3VY1MwKZN0ObPZTWghS1YFGPlQjjtAnsu1ZD/pINfl9HjItDGJTA+Pcq+KyqkPfG/4Nz7r+dlH9s8D27VwC8goXisShrt/p1jv6Ok07jMCeBDVV4w+/gK36WA7COI+nqdBqub/aU9ZSbvqKyu8pZzwxRG3ZDU2BqzRoaGzFU17CiNpepqZUsf/8pwhoV/1ICBmBAdCY0QvnZc8XNuq6MvxjQOMDv6ln+ouGKUhJfVoPSWs6faikO7TOcjLYHf11G1p+7VxH8TfJr9o8ImE1CEBATEjA0NIRUUJtL9PGz4w1OXFJC/OveZTs76TQOc/r5gKAXOt4bfmahAaBtcVDJy1s4GQyGPaJOK8zc/E3yazh+go6zisGUJNldigC3MHmKLrhzVugJDauVJJu+VYw++4M4gzu/FGv0jLdy9dx5KzRcod9qHAOrKj0v6PmcPpDcOudOdli02R844qdNS8VQWwc4djbrbVhrQD35BXoVyBubl5g0yKbQlKO288fuddJpHKdd5eN8TSN6VYvTV3jCTzssi0N3lLmsCu4pxEGJHu/rc9s5QKCFBjjmt3dpDjHfJBD1tef3GUhYa0A9+QVakPUUGq48jntWnHfUdkJYYGxRNpwCfsYgYf1y36zdVXdOtdSXEBPiA0nJBp7wMx47Ttb8CsSMNL+v1/MBcgVf26630HPZ0xG/oZ+1UdcW4/25AyCUeyLY7eeNYHLErTeiwPut4CiY7n1xX2HSaN6/7kHL3N7YIwmwNi1wfv6e8FM6O2l8yr/s5L7Al7brTez5axF7H+0eKR3x0/ywiYizDtFxkXejZ29Ma4Ldft7cQ7C49VvBUb3dh8AvQeC6S26wfLXuEDEvh0Pn5VoqxDtz/gkkP6Wjg5jlnpbeDRxccrOaxtUsK6P6zrJe51P41wqGWyXnccZP6ep0GCFsSTswLIuaZb3P16dnL0g43O7fc+sp+q3gSEjzzIJvzn4Faq5LzQ+bHO538KFIhDPqacpQbSC7/z4azQT7mAltRnpA+fWFoc8ZN010tKWMJIDYoTD13C29zkd/uMbGe9PTtrMcnxCp/v+HgjGsd6eg4D2/YCKB4CQz6reCo6vDsxdOyHLzoptG2PY98ei/TSLxrEPs/vskMt/Voo9VpbeN+7OHhiZP+fUFnHEztrTYCLKkbe10GIOfssXbthNWbAZBwPCnFLJe3MGOh/N7iZmKUO7bjuPBycfRbxP5dDS7X17yqHaJ6bfcO1Yi6MJo3zGekfsbMOwp7/bvsIqk7Rw6CE35ARAEtKlDnCbN9YSfp+i55OYvPOWm+X4jdTODo/paw6e2UxTEbzdiAEb+ag+96ZLV0e65r01vB+z1RCCfO1fotxpH0lD3y3KGbbtg6ljL9KLxItUA56zymNLVSfgna9XAL2s3bitHKXOB5bprS9h3vfMweU/4eYqukX66pPeAN9yCXQkeXPNz6WSmKKoTmZNK7oFCUobnU5VgT0UD+dy5Qr8VHBUbPavhweqtlunIvLv+R+VbY9Gu3YVmXKFf14855HpM85ifB3Bml/EVfnPTiDZ2BDFpUECXPV3xM7b1bSIkgIoNwV0F8waBfO5cod8KjjGzRrj8XRMRgTY3GxQFfVU1AKtnJDMzew/GcXlUzR3kl6twxAdrGPa7lU5/d8cvIPAioY41zNwErdancxjLxmJo7M4ydeL1BLWwdoDgqu0CFfnrD4LStz4iWNz6rcu5q9qxmuhoNAnxtBekEbZiG8fPnUDcf1Q/ffMyqyY8nM5JeWjaDegqa22W+cT8XAx7yxFHjUApP4ixrc3raMtg1Bf1JBBMKB6LUaexqZ1izU0cXYA+PsLj2iqOoE1LRYmJsuQM9Re92Xad84oJ+9S/lIihVjvWGsGqHeuXfilJUgLwIjAGUICrgF3Am0A2UAFIsiw3mPa/G7gaMAA3ybL8mWn7JOBlIBL4GFgmy7LLN3X2dWV88vh3Dn9L/lLk9MR1PPzy+WjHTCTt22OWzOHmObuho4Ow9XsxNDWhpCTbHK9U16DNSKdhbCKJza0YK92ofw4iJF3xCxQ8CQQrfG4Hy9dOZIRVas3Z15Xx+WeHMOzYg2HbLrRJg/wyJuoP19jWnfUTvdl2UZsP4q/VwRE/v3OqBAjBeO7A/6nKE8CnsiwXAuOBHcBdwFeyLOcDX5m+I0nSKGAxMBqYBzwjSZJZT34WuBbIN/3Nc3dhV43z6rD/cXZ0OZtvfIoVdzxO67A4xER7o5GhsdGmBogZxtZW9FXVxL6xiq6MQfbahlm9d5H0OBidZwPr+Burz3vOS2fE0m4ns855xXzytx/Yc+Vgy7aeKzZHbipDOywLQRfmsSOcs4xUnqBntq/ebLtAlI50xC8UhAYE77nzWXBIkhQHnAq8BCDLcqcsy8eBhcArpt1eAc4xfV4IvCHLcocsy+XAXmCKJElpQJwsyytNWsarVsc4xcJ75rj8fea6n3PNwRmMWX4jUZVNHL5oJMeuKkVMTETQatEOy0Kb2l3gyNq4p4mOtvzXbq+wbBdHmvwDeggKbc4wmi+YamMwdMfPEfSnT3L5uyPhZ4a145Z2aGb3OSsqbfgeuNjAwrtnkbbCuY7RlqKgP3gIpavT71WV1kVTvQ4qPOe3bseNPoUvfRssBIubPxpHLlAL/FOSpI2SJL0oSVI0MESW5cMApv/mQhwZgLU5usq0LcP0ued2O0iSdK0kSeskSVq3bds2ho5LJ78sm1Ez80gvTKF40TiO/Wcq9229h1MzyrmAMcTtErn1+XhOTOngsvFD6ZwynBm/O5vMfxuYeN440gtTGDUzj7S7ZpN062xy/raA+FNHM31JMeFpg5n12zMQtFrO+d08hLYOS8fMWTqNqFgdpYuLSE4Op6zxBMPmjyV3chbj5hay4j8bKF1cRFR8BHOWTkMTEWE51vx//s0z0EXomL6kmMSMeIqj2h3eU1xyDDOvKUUTG23Jt7Dw0QsANf+CoBE47Ypi4pJjKF40jpTIDkbNzCO/LJuh49IpOns0iRnxTF9SzOil2+ls6SDmh312fOYsnUZUfAQX/dRO8rAExs0tJHdyluWeknOSbO7J+lhn91QyfBBDx6bZ3VPHHTOY8tz5iBNGcvat0xF0YZbsVYbOLgSNwMxrSi33ZO4nR/eki9Ax/+YZDvmY76l0cRHJOUl+31NSfhqVWw657CdBI3T3k+mezP3U2/dUvv6gT/1UdPZoh/fkDD4bRyVJmgysAk6RZXm1JElPAI3AjbIsJ1jt1yDLcqIkSU8DK2VZfs20/SVUe0YlcL8sy7NN26cDd8iyvMDV9Z/6zQvK8r98Ybd99z8nMeRLHQk7m/jtW//i86ax/D55GwbFyLydC6lpjGVm1h4+/HESSZsEoo/qifhsI+X3FrPqikdIFKM4ZcsiYuaphr7MVTH89MxYjo2BQT9B0sojHhkBF94zBws/QUAcno3Q3Iq+Rq12ps1I97kyG/g2pzZXhbfmZl0pvreuaw1zbojr9+zlnOhmPmqN4KFll9mUY7BpuwBA0GrRZGdh2FuOZlyhpYaMOzhrm4X3zGH5/V8ihIX1iZ+LKwS67XojH0cVUCXLsjnP3dvAROCIafqB6f9Rq/2ti4dmAodM2zMdbHcJZ40ze9ROBCNM/ccmSsLh98lqkZ9mpYPlBe+Q9pCOFn04aT8o1E3v4sPnn0T7ZTJdiUYeqithd1cL0b+P7b7Jm3IwSvVcOvc76ue2o0R5Nue34acoGPaWW4QG4DD7uDfw5eU1L5lac6tc6rg+hyMIxWMRoiK9vq7NOUQN2swM5kepNWHPimpnzH228TCBfPDB5ITVpEZCN4xNQCzI8+g4ZwJ1+V++YPezxWhiou2vdcoEn3kGAoFuO2fwWXDIslwDHJQkyVzSfBawHXgfuNy07XJguenz+8BiSZLCJUnKQTWCrjFNZ5okSSqRJEkAllgd4xRmFcwa4pAUqk6HQf+rZPlLMxAFDXWGFnI+uYYF2y5m8nM3IzZ38ETml/znkUfYOfdZDCjsXplN2new8dR4brzgOoSVVkuTa7YyaMEeVk8fjLg/AuPmHWizMjm+pNSh4VAcnOSUnzX6YqQy52Ww4TbFcy9Lza4DdoZkrzm0t6OvquaBuvEAXFpxGnuvz1edykxeoe7azrcLq5p1grwBod2/tp+zdBoj761wGAagPeZ/7gsxaRD1V5f6dGyvtJ0D+OXHIUnSBNTl2DBgP3AlqjCSgaGo05ALZFk+Ztr/16hLtnrgZlmWPzFtn0z3cuwnqNMdl8QWJl6mtJ6w9SKc/VMTX46J5ciNZYgdCm1DBIas6STss3VqFqWCHIybd/DLPbtpMkby21ULSftQR8xbqzHMnMixgnCSn3Ps1NV4cQlH5nUSsymC5mFGZpVspc2go+7UZnUkVxTQiHSeMRFFgMQVW+nJL1QQFR9h4eZvOQVfUwceX1JKa4pAxhPr1OtrRKrvmErGAyts+AUMASiXYIYzfuLIfBomJFl8hlyiZBysch55LI4ucFpkyxduvsLZVKXfOoD97uK/KCvf2Giz7cTHecSfuVcNCjvWoD4ogoDx1AlkP7ib0xJ28ttVCwnfH0GOXItyoBqlsxNFr+f6PXt5dkQ+CBrEkXl2naaJjsbY0sLu56YQltjO4sL1/D55G6NXXkJsZAdRj8ajbdUjrNiMNnUIxaels2pPp6UeaiihdHERPdtOm5WJ/qB37sratFQOXpRL2qO29VcP31ZG1r/3q1MzL54vc9Z0R/yc7dsXcMav8eISYivaHDrTCUWj4ac9FiHdW1ncPWk7d7Cu/nfS5Rzdu7rSblvi+appxFB/zPLAVv62lLa71QdscUwt++f8g/Y0PYYde9AkJ6FMHoVm/EgeueMSmDqWlkWTue/D16i/ptSmUplZzR9x3RpyluxiTVkCed9cybNFryOPeZmOQTq0DeqDoK85wt7VlYjH+r4wsTY3226bo7bzVmgAdA5PtRMaYnIyW257hqu+W0HzBVO9Op9ZEDjid/A33Ql6BF0YM7733bDsLxzxA4j79yqnHri7roui6YNuU15vZXF3xs0bKM3up1v9VnCkjUi22+aoM3Ke3UPMwmoOnqpn9vZzAShYpnbu7V9/xF/+8yJPvP8iUe+tpqYkhqOTNPx2/7kYFjSQf/d2u/OBap8QdFqGX7KRP5Yv4KA+ilcfepidSwfRtLjEwk9ffsApf5cp8M1+DwFIqqzfX2G3zVHb+QJHwXeHLspn2pZF3Ln8EmLesq8Po81WI32FYufFkXvyE8LDSdxttBwvFOTy0aExiENSHB3e6/Ck/XpG8Sat03KofLCTvQnYvQSib43t7qc6/VZwNNd7ZoQy1NaidHSgdHQQ8asodne1UHv5RA7fWsZpkUYmhYcxTKt6SKY+voK8B3ey93Ayabd38fV25xG07VNUZzDNnEP8ZfYirvnFLUwu2otxSR2aiAi3/Az5mXbbxBHD0WZmsP+vJQiTx1D+5xK02UMtNVUDBU/bzhqaiAiHEbBK6Xib70P+toLo+eUMv32lej/Dsmx+X/blJ2giIqic271yJSYn2xiae/JTuvTEvK0u1xqqDrH7qgSibwvvNtQGIWu9NTxpv6bZI20c9oZ8XsWI69c43FeblYnRbGj181586Vtf0G8Fhy8wbt7BdVffRFMObPnVMwAYFCOn33wDmkg1Bb0xO53YNZHsvSKFvGFHbE9g1am6z9V6qmLuUPT7Kwj7dC1NZ7QR8VQi+T8Y3Ltq9zCMiQnxGHbvQ19VTf7L9dROjCX2AByZk0Hd4vGQmmzJrelrVKw/OHJFEbuemmizreOsYqIecODCbZomCo3NKDpbYXNGVBeH3sjh2os+tmzTj8hQV5mcvTRGg8X7VdHryX2vA+O2XSinTFBTQwq2j7FZq+lLxH6/H5TuhD+uElIb649Z8nbsecK/Qt7BEqL9NgNYTJL9Gron0H25HuOVRTx6LJdhYXXc9d4l6AoE4j8zLdc9dZhPc15nX1czw3UxnBl3qt16vjg4icoXUokK76RpbTJD71PjNIytrYR/vJb/5ZUx/fRE6kclk/LmNo8crKwNfYbtuxlsPUsSBAyKgjbT5FBr7fLuQwlCX9rurbsfYsHa62y2hR3vov3WwcARC09rY6ijuJD1HZ1ckLuR5+UzGYpqH7HkGzUd646f5nvV+GfUaRCspqfmFaKu9ESEispeK8/oSfsZams9Pp/1FLvw93sgMRFDQ4NPK16Gi8bB2speL6HZb1dVLh5+nVJb7rsTlRgXhxAXy9HnollV9Aa31UxhxyQ92tQhPLvmHc5afy1nDN3J52+WkP6grQFQExXF3nvH8+nih5j7w43kX7cXISoSpbXN4pg1pDCNIwnJiHurMdTVI+bnIhiMDm0OgYKg1XYvDbtAck4S3radUDwWZf12v19EsSAP5dARlw5svvADaL5gqkO7SqDhKz+P4cfS8ZBRmRzZXkX8D0mcmOY/x5NuVSVvqn/qqKGxEX1VNQkPRDHhiRvYUJcFgoC+5gjXnnc9WZdXsf3USKJqFLus2cbWVnLvXMkZ7/2K9Ld0aBLi2flQJvv+nmOxA+ROSIVVWyzLWoa95SgNqlahiY52WZ0LVHuHtzCru2KyawOZL22nrN0akNHbsGuvW6/XnvzM0bO115dydGkZHfOLHR7XFaWxq6UaKFjX3PH32XMLPwbz3HHJGKdNYN1mz7xjfUW/1TgcOYD5A+uksuYYBXfxJB1nFvPI009zT84Uqu4uY+ai9ewtNaDo9aojTmOH/UMgCOrIKNs6CYn5uTYxMNqcYS5XZVzCjYreKw5WPWE9ano5ZejJTxMVRcOicbz+54fJ1IZTa+jg56PnWwSQmJysTg16sXL8gT+UWjK+BaX9fIRbbj20GXf+MCedxnHKJZOBwJX4swiNpEGqTUIQOPXTPS6PCf94LTffeiMAmfevYP+cCMt5Trlksp3Q0ERFIWh1dkID7PNZ9BQaXt2nm5fnlEsmexyv4QvEUSPQpneXtVSmeh4PA919a4EgEFFvoNYQyZjXb+LS626x0Vos9oReEhoAOX/tNmbb8QshuOVm9UwKWi1VL6dz/DLv3dv7rcYxR3OBAqDNzLDkFA0kxIR4Kq8bTcYDK9zvjLpc6Wr9+9CvylBEGPrSLvsANw/mtMeuKmXQP5znOPUWPqXt93Du3TG/mPBPuqNdGy8q8cwN2x2+ykSz4FivGv3+v+HEpSWEnzAS8YHjpeKTTuMw5xVomuQwdYffUNo7GPr2YbtlVc2EUQ59Kg7e3L1UKejC7BKq/HDzI6y58XEa5uTbzcM1kZEwZSyCLsypt2XS6xt8vRU7LLxnjm9p+wXPHhdroQF4LTScJqOZczgkhMbJlMgn/rVVToWGK/R7jaM3oT99EmGrdtg+rBqRPS+PJ3+J7Yts7d+f+OMgGi+Pt5l+dMwvpqFQR9wBA2K7kYjPN9q+vKb5eW9pUH2Fnvcjji5A0Wowbt4RXCKCQNvPiolc7v1L8v8ZJ63G0ZvQfr3efoQzGsh/xn60tp5+RIpdnPGsbdGn8E/WkvrYCpJv3o8hXEPHJxnsfbSEzrmTLecFfBcaXjiFBXPEtL4fTVQUhm27qDwr0S7PqDV84SdMGu3yd+2wLE5c6XkOk0DzcwkP+86T/K/B6tsBjaMnPJjHa7OHqrk8HaDlvKl0RWkY/E2l5aWxdsZpurAEow7+cN9LXPfpVWS/b1C9UP2EUjqevRdHkH9j7/sx+Apt6hD0NUcQE+JpnlFAc5pI8qYWhPU73To6tZ89hYgPHWsLYlwcS9Zu5Z8FwyzbhPBwu5wn+tMnof16fc/DLehtpyln8NTeFOhSoPZE7J/9k07jMOdjDDg8EKTOhAZAwwiRhNdWMef87lWLnY93j4axb64i4Y11PP6zcxlxy0afhIYmNtZum7ByMxnfeHa8y7bTiP67tDtxe9bXHEEzYRT7bx3FsZFamofCLa+9ifaLJEt5hfprSpn+6oWAbdX5tusbnF6uc2Iez916vs02R4mSXAkN8DxiNdDPnqf2Jk+Ehl/cvFAi+q3G4aogk88IoB+AdWEcR6OfP7D4LfSEh6seror2aHOz2X7PYEY+dBxjxUGfeIsml2kLLWvXaUFAzMvBkBiNZkcFlTeOpa2wncLbDnL4gnwWX/8FP1xbSueqrcR/P4gTt6bDmq0BTcTjL8Liozlw9Xi7lAKhgGAVZOq3GkeJNCHwJ/VBaGhzhjkcYa35mV8+Zx6PzmBOQ9gTTuMgPHyxXLWdfn8F8ZvDOPizFAQfA6ashQZgCSAE1Pyre/ajrTyKsamJzL+soOC6nRhqaxny0nq+PFrI1PxwapaVUt8ezdFik3bl4YqOGbXX+5Z6zx3E5GRKL5xI5me96HLuB3rlvXCAfis4fvpqd19TAMDwkp6Dv7Z/SHvyE5MGUXlhd7SkEB7eHbTm7Nx+JjR2BhtuDoTDkKdXk/7QCo/yMngCR0F+1ombzVMEpaMDzZxD/PTVbgZv7SD8knZSnlmpphUo9M4F/0SBGyHqq1CsP8aWj7YgtLRZtlmW163PqRHVxMUBjFZ1F6YAwXsv+q3gyC6yz2fRF+j6UyrDHrPP+tSTX8X1heRf0T3H3vXcGG789ste5+cINtwcaSm96IHpFkYD2UWZaL9ar+bbUBSM7R0oFd5lKMu72Y3viJcaTE9+1mHy5tUOMT/Xsq1m2VT2XBaGxsNKeJ7AUFTgdp9gvRf9VnDUVzo3lgUEHo4U2q/XW9IKWqMnv6w/2c6Hs/4r8njl7O7zZDnp8F7Ir+Cs7QLlvu8vevI79KupGNs7bALN/IYfwrG+ssFG4JqnZobd+7pPr1PTTNpobX4ancUtexHCwzFOm+CaWxDQbwVHeIzzdXZvIMbFOYyoFGNjXa7lu4M7fjHf7UKZrearELRaujId2zM8sVtoM9IDws1RnZCgQyPa8esoagHFiJIYuCxoPdHTnuSqTT159qxTMZgFspiX7Rs5E4wtLYgpyQ5TNnrDrSd8SVvYbwWHLjwwo6OxtRVje4ddKLqhsdGvsgHu+BmOn7DJaqWr9DzxixnaHNVvwRzBq4mKQpg8Rl25sMpp2nNu7IxbX2UNt4aYGG/HL+eSn1Sj6nYH8/cAZUPrqW25si95++xpctQwfGuNxFfYJZXuoZHqwrVea6nKCc8r+ZnRbwXH8cNWN+uHOq/JHQZGg1cZmzyBDT8P4Es5yJ4RtPUXjOd4YQwUj6H8OdXw2nRhCXv+NEHNdi4ICFqt19wCgT1PT1VribiBof6YPT8X04rGj7L9ZKbC2lgLrgtmedt+DkuGCoJdnhef0EMjPX640aYAeU840qR8MYL3W8GRMcpqvuvl+r42K1PtuBHDnY4C/iaEseHXi7BOiJv4ykri/70W1mxl2KXq6Jy45jA5/20HRUHMy6Ft/kSvuGliYx06nHmLgr83IRg86ydv+MX82X9u3iIQfVt/dQmKweh+Ry+RMSrVzqHM2h1/32NJ1F7n/1J1vxUcu36s8PlY5UQjglaHUnXY8Q4a0cYY54vRcNf37gtTBwKGgiz2PWT1IBgNCJNGW0ZMffkBND9sQl9+AMPeciI+WOMVN2HIYL+KTFtobdquZhHzAN7wczXfN0NMGuTc+OwDAtG3Q76sDki79oQjbgfnxls0jdVlL9AwvluD89Ug3m8Fx6SfuQ5qMsORgdNsv3DqYmzKqq3ERqvqfaT3hZYnLXReNySQEDbuIv+14zYBUJo9jjNqiybNwRtuPRMMBQOBbjtjcws7/jgkYOcLBD9XYQv+wBG3rAfXYGw4DkDps7cx6i/qtLjtnCmU3+edU6IZ/dbl/AytpCimQsK9ntHZh6Q3gkbAzC/UEMrcoHf4OXXT9wGh3H7ecBMT4jE2t6Do9RinF1myx1vjpHM5P/v20y2flU7fVz9cQZuqjlK+JL2x5hdMeGKb8Yibn/4j/kwNeqPtAmn87qu+9QTecDMcP2F5tivPiPDIM9WMfqtx9HZYvaDVqpLbaAipACtvIBSNRtnoXdFrs3allI5HWOm4Dqo76E+fRHjV8YAsPw4gONCMH4mmsdWyUld1dxmZ9684+TSOhXfPdr+TL1AUOuYXq5LYtAyoiYnx2ojUa/y8gLC7wuH2hXfPpvb6UofaiXkE8lVoAEx9ZK1fQsOntjNpSN700+5/qEmUnAUTOkMo9K0z+MJNjIvDuHW3zfL+sOdcZ2gb0DgcQBMba2vx1oho04b45GsRiqi8rwxNB0QdUQKaANmMnmH1vQnztczlJLRpqQ4ryDmCftYktF+ttyQYGoA9nGkcfrlfSpJ0C3ANoABbgSuBKOBNIBuoACRZlhtM+98NXA0YgJtkWf7MtH0S8DIQCXwMLJNl2aVgWHDnLD7461f+0HeKnstkgih6LTR6k5+/WBqfySd/+hLjiUZ6Q/r6KzS8aTvztRRRVZ6VxDiE2joUgwFxeDadmYmI326wCAlrmL97KzRCuW895ubn9NvnqYokSRnATcBkWZbHACKwGLgL+EqW5XzgK9N3JEkaZfp9NDAPeEaSJLOL27PAtUC+6W+eu+t/+NDXvlL3CNZqvC+u573Nzx98dNtbGKwKHYcafGk7w74K9f/23Rz5xRS06WkY9pbTFRP4wL1Q7ltPubnyLvUE/to4tECkJElaVE3jELAQeMX0+yvAOabPC4E3ZFnukGW5HNgLTJEkKQ2Ik2V5pUnLeNXqGKc47aoSP6m7hr/Lu73Nzx94y233c35WULeCJtp9IJ1PbWc1eqY8vcKiIUYdVLXH8M0V3p/TCU6GvvV30PBZcMiyXA08DFQCh4ETsix/DgyRZfmwaZ/DgDn0LgOw9kyqMm3LMH3uud0OkiRdK0nSOkmS1tU0VjN0XDr5ZdmMmplHemEKxYvGEZccw8xrShE0AgvunAV0G4wW3DkLQSMw85pS4pJjKF40jvTCFEbNzCO/LJuh49IpOns0iRnxTF9SjC5CZ8nhaM4ebf4/Z+k0ouIjKF1cRHJOEuPmFpI7OYvcyVmMm1tIxYaDlC4uIio+gjlLpzk8x/ybZ6CL0DF9STGJGfEUnT06KPcUmxzt1T2V1CQwbm4hQwrT/L6ntJwEt/eUmBEfsH6adUoSEXOLmXrRRFLHDmXCpdNs+ik5J8nre2pv7gjZZ6/leKtP9+Ts2XMGn42jkiQlAu8AFwLHgbeAt4GnZFlOsNqvQZblREmSngZWyrL8mmn7S6j2jErgflmWZ5u2TwfukGV5gavr33P+H5W1725xtUufonjROEKVn6/cjDOK0HzXw0moF+q1BrztNCIoRrQ5wzh4bjrp/2vy2P09KPwCiEBz643l2NlAuSzLtbIsdwHvAmXAEdP0A9P/o6b9q4Asq+MzUac2VabPPbe7RPV2zyznfYW+5ieOGuH0N1+5mYVGx1mqm7IYF4eY6HuEp5gQjyYiwnI+f/k5hdEAioJ+fwVpj6xAs0d19/Y1eK+v+9YVAsrNRcoCfwRHJVAiSVKUJEkCMAvYAbwPXG7a53Jguenz+8BiSZLCJUnKQTWCrjFNZ5okSSoxnWeJ1TFOkZDWe0ldAgHlWjcxAAHKI+EULjw/vW07c24Pi4+E0bRk3daOof6Y22JIzrDvV6MQP03k8aee8ouftzDnHVEKsz0+xjrmKWSfPUEILDcXmqQ/No7VqFOTDahLsRrgBeABYI4kSXuAOabvyLK8DZCB7cCnwC9lWTYzux54EdVgug/4xN31uzpCc0XADOGYG/W9l/N6Grbtcvqbt21nruehTB4FqFXpBK0WNKpw0uz3zb8l7+lyIsQuLnnhFr/4+QqvpiuabkEcss+eogSNW791ALt6wjKlckvoOmQNHZdOqPILJW42NVdMcMdPmzMM45Fay8pX40UlXhe29geh1H49EWhuJ53LedLQRPc79SFCmV8ocXPkI+OOX/lDsTbL5Ylf7Ak4L1cIpfbriWBx67eCo2Kjd+nyg53B21t+wUQoc4Me/HrYagStluyrKm0jOZ2UOjDOKHK43V+EcvsFi1u/FRxjZjlfNXCEYHtJOuIXiBR8gUBPbgEtOxAAjJ1TaPks9mgzRa/H0NhoyZO57+ESsj5qtks2DSCu2t4r/Lx99oKJYHHrtzaOXqkd6w5e+PcHuoanJ/A0oZEdt17wxfAVYlwcWp1AR71nGdeN0yawTwon/6bV9ucqyMOwa2+gKfZJ33qKgdqxbjD7ujLHP5iKGvfcFgiIgwd7vK9Tfr0IIdKzRCx23EJEaICa1vH0yyZ4vL/mh00OhQaAcrB3DJh90beeIljc+q3GcUbYRYomKgpE0SYaUxMRYZmWWKYnITSiDmAA/Qknncbxs7+eizB4EEKEVZLe6GjaZo6FopEgit1FiYIgNIRi2ySx5piAUEQocdPmZtttCyV+jhDK/ILFrd9qHKeffr8SvuuQTdIW7bAsak/LpCNRIOOjGseFcHoJ3iSQGUA3Kv5USvZvAp9MaACBwUmnccz9RaHlRdVERIAgYDh0hPh9bUTWGYOe1r+n0BgYlTxD+vf2hjxf+flbRMtThFL79cSAxuEG8wruVJTqGoxtbZTfX8KQNUZit9WjHxyDIVxE+/V6h16JAzh50dtlMv4/4qTTOE6fk47SpQdFIfKwQERtJxgMGCJEwlbvBEDQBdfpyxrmPAihiFDmBr7zC5bQCOX2Cxa3fqtxLEy8TGk9oToBnbikhIQ316l5JvNzQyItf1R8BGZ+oYZQ5gYD/PxBoLmddBrH+PkjLZ/jX1+lLr0qSkgIDbDl5wmC6VXqLbdg46Tn14spFYLVdv1WcOxd3Tu1N63haKnQU3jLT2lrA40YlJiaYLSdPzjZ+WmsXAgCjWC1Xb8VHGkj7GMTAg0l0vsO7pxXjHZYFmlzvQuwUjUmY1BiaoLRdv7gZOfXm7aYYLVdvxUczfUtvX4NV8lwzBC0WvXPVC1eEWDP9ZnsmuqDm3uQ7E3BaDt/MMDPdwSLW98tO5wkMJSORdvUgXGTGokZ8cVGcj7Rk1E8FJ/dz/pprdoB/P9Bv9U4YpLc1+cIBnRb91uEBnTHx8QM8s0ZSRwxHDFpkK2twzpILwABe73Rdp7US/EUodK3zuANPzEuLqi5YILVdv1W4zi8uzYwJ/JydBfj4lD0eoSoSDXA7oiaxF0zfiSauhMoUREoNbUcPeTBPNZR8F19A4aGE7bbe/AT83P9cqc3t12oOkwFrG97Cd7wMzQ1+V01zRsEq+36rcaRN3Wo9weZRmtxZD5oRDQREa5D5TUiQni4muFaENAOy2Lv3aMpv2s8ldcU0HRKDpqICARdGO1pMRgHx9NSOBhNfBy5xVmW65mzVWlzsxF0YZbVGuMp4+w0CEP9MddBeYqCUu1fTIy57YQwnV/nsYaxJXBz67zSHPc79SG8evYUJahJpHx6L3xAvxUcmz/Z4fUxFsl/TB3RjR0dHD0nz+m6ujZlMJrwcISRuTQsKaEzezBpKwx0JhlpGdFJeEMXaDSIKYPRNXbCrnLC6zo4XpbF1nVHYepYxCEpYLrukVlptM8ejxIZjiY2lupTIzn4m1Kv78NfLcHcduYyAaEGX/o2mAhlfsHi1m8FxymXTPZqf01UlEXym6cXKApJL60Bo4HWc6eizcq0ESKGhuMYGhvRHGsi8dVVaDoMRH6+mcL79nPW2K3oo0SE9CHsXjYMbUMrwrBMxJ/205KqYexNE9AebeToguEcunY8lb8r43iBgqBA199aafjZaBQRcl6uRJupVrw0zJyIoNWqwiaAsK4JAt63XbBxykW9kys0UAjl9gsWt34rOL545gev9jcvl5qhTR0CgMakrrekalCaWxBE0RJlqXR0ANCRN0S1M6zeitLRgfH4Cbb8eTzR6ys5eloqQz/tpKEoCY7UYhyZzbTL1vPfDdWcunwbLRkCQ9a0YhjVzJ8WvEnNVB37t2bQGSuQ+89KWkenUT8jC23OMPafp3KpujSvB3n/DKJiVrrNd2/bLtgY4Oc7gsWt3woOr8KHBcGSJUxMTETQalHiYgDYd686uqWsaWTXvQVqNK3ZNmFyA9eu3Kaex2SkNE4ZRdSBFvQ1Rzg+UiG86jjHRgsYGpthzVbWPjmRq8YO5XBnPGf+bBWznvuRQR9G8dcnLiJ1dReCHlK/qUVpaiLyYCMxhzopvzSD7PcNKHo9aY+ssOXv59Ksfn+FzfdQDguHAX7+YCCs3g3maC5wTdxqtUTQam0MVOLgJJTWNoytrQjh4SgdHXTOKyZq11H05QcQExNVQeMk5eCh28tIf3glbT8rJvrLbaDRYGxqsvzect5U4tdUU33OUITZxzhxIJ6UNQKJW45j3LoLQavrFlAmjprYWIQwnaVq2v8naCIiLFnLAwUxLweO1mNobBxIr+AHTrogN3eS1XoJrKdV21BXbzEwmqcjusZO9OUHADhygSk9v5PVjfSHVoCiEP31DowtLTZCAyD6ndWcddlIhjy1kiHn7yd/2RriX1uFcYsa7m95iK2EtrGpyXuh4eMUJtRGzBMLJ9h8DwQ/48FDGE19G2ihEWrtZ40BjcMN3GocjmAa4bVpqehrjjifAmhEBI3QK8toPbWfAfjvlzKA3sNJp3HMv3mG9wcpCoIuDENtHS2LpjgfsY0Gv19uZ/yCITR6rqL0hE9t14voKTRCjV9P9CW/nkb+nggWt36rcfhbkEmMi6Pm4tEkP79K1TwCHB9yshXtEQcnYairDxgHbeoQVetzAI/49WHJi77sW4caq9WzO1CQyQ1KpAl+HW9obGTIq5tpP7sYTVQU2vQ023qkJjgqLdgT2ox0u23+8utN+MLNcOx4QDk4ExoQ2m0HfcvPocZqNeAFi5tbjUOSpH8AZwNHZVkeY9o2CHgTyAYqAEmW5QbTb3cDVwMG4CZZlj8zbZ8EvAxEAh8Dy2RZViRJCgdeBSYB9cCFsixXuCMuZV2jNFT75vloLbXNqyoIAprISBSDwWIwBTwa2TTR0XYu14kZ8fjKL9DoGZMSStwcYYCf7wg0N380jpeBeT223QV8JctyPvCV6TuSJI0CFgOjTcc8I0mSeXnjWeBaIN/0Zz7n1UCDLMt5wGPAXz25oeyiTE92c4yxBZaPFiGhKCCKtkIDXAoN83zTUZyGX/w8gYcrKpqoKJrm2xaL6nVufmKAn+8IFje3gkOW5f8BPdcJFwKvmD6/Apxjtf0NWZY7ZFkuB/YCUyRJSgPiZFleKcuygqphnOPgXG8DsyRJcvtW1Fc2uNvFBuLgJED1GdCUV1m2Wxubei6ruoOdkPGDnyewCc/21B5jNBL93jqbTb3BLZAIBD93BmJ1J9+Ws0O5/YLFzVcbxxBZlg8DmP6bgysygINW+1WZtmWYPvfcbnOMLMt64ASQ5OiikiRdK0nSOkmS1iVNiWbouHTyy7IZNTOP9MIUiheNIy45hpnXlCJoBBbcOQuAhXfPxlB/jAV3zkLp7ODU80cRlxxD8aJxpOXEM2p2Ifll2Qwdl07R2aNJzIhn+pJidJFhFiu1eX3c/H/O0mlExUdQuriI5Jwkxs0tJHdyFrmTsxg3t5CU4UmULi4iKj7CkrK+5znm3zwDXYSO6UuKScyIp+js0S7vSRMdZXNPAAvunIWgEZh5TanlntILUxg1M4/8smwyRwyi6MzC7nuK0FF28USf7ik5J/D31LOfAKZfXuzynpJvnGXppxk/n4YuQseZvz3bhsfPbp/h9p5S/jXfp3saOTPP63ty1092z16Ezqdnr2B6bkD7yRk8WlWRJCkb+NDKxnFcluUEq98bZFlOlCTpaWClLMuvmba/hGrPqATul2V5tmn7dOAOWZYXSJK0DZgry3KV6bd9wBRZll2a8JfNulvZ/s1em22d84oJ+3St2/uxg/WKiiCAoEETGYHS2eWz89ComXn05OeSQhC9G73lFmy443f41jLSHlXd8rU5wyyOe95CHJLSHfAYQH59iUBzC/SqyhHT9APTf3PrVwFZVvtlAodM2zMdbLc5RpIkLRCP/dTIDscPN9ptC/tsnYM9VYijRgCgzR5K57xi23IEVsJTEx6uZuCKiYbxI9zRsKCnamzHzyrqtlkqsTu26ZzgRYQ6artQgjt+aY9115pV3OQU0cTGqlm4HExdfBEanvDrSwSLm6+C433gctPny4HlVtsXS5IULklSDqoRdI1pOtMkSVKJyX6xpMcx5nOdD3xtsoO4RMaoVPuNDrQnc6Sr0K6O5obDRzhWqKO9rMBuXwBNeiqG2lqUIYMQqzzPptRTW7DhJwg0n19M7XWlaFOH0JCvNrtxRpGq4YwfQfy3wfOcdNh2IQS3/Kz62bDL8ehqFhT1541B+0EUdZdPCjg/7bAsxIR4u2v2JYLVt24FhyRJ/wFWAgWSJFVJknQ18AAwR5KkPcAc03dkWd4GyMB24FPgl7Ism5clrgdeRDWY7gM+MW1/CUiSJGkvcCumFRp32PW9Zy+aeRlSv78CbfZQBEEg9fEVTrUT4xFVWBi37HTpa+AIYmIi2mFZdvy0qUM4Vqjh5pvfomFGDsM+NEXqrtqOOGoEyrqfwIWhNdDwtO36CoHgp3R1Ihbk0RkrEKPtoDFX3a5N8//FMvPbd3UmB/+ZYdmu+dxFNrkgIVh92289R/9y7aPKNy+udL+jFWxiInohk7g4JAVDbT0YDcy8phQzv/YFU6gfpeXbGx5iyvu3EnVQRNcIKc+scHPG3oE1t1BEoPiJCfFgVBDi49AfrELQhfHzbTv5x4wy9Idr0OZmW1IO2EToakSOLp1KylOO+8ean5g0CEP9MbRZmTSPTyfiwzV+8/YHge5bZzaOfis4ztBKimL0n7u54wMNQSOgGBWLg1nHWcV0xoh0XXaM5FsNtD2tR/e7BISVm9n992JG/NwHo66f3EIVvcXv8K1lrL/tSRaOnImh0dYWYO0CX/3uaB4e9zZLP76C/JtWYzhtIuK3G1zzC9BAZHFI9PX4ALfdSedyfvbtp3u+s0Z0umZvbO523tKmDoES50tQnqJZKmHBPXMRkwZZHoLwj9YS9846ki+rxbBnP41vpzPl2Q3seWoq5Wf9HTE/1+/regqv2q4P0Fv8Mv+xjel3/NJOaICtC/z6qS9Tb4jh32c/A4JA2CFbT0yH/AI1ABv8i78JVt/2W43Dk7B6bVoqxmb7fBkWlIyDNdss3qGCVkvzOZOIfnu1+ruPo0izVIKuyUBripbEV1c5PUfNsjJahhrZsfhpPmqN59n8PIf7DcAxtBnp6KsPud/RS9RfXUpTLmCEtFV6Ir/d7nUWd01srNcOhaGIk07jMDvWuIL+SC3Gllano7m484CNS7mi13cLDfB5FImRVyFNiCdpfYNNQiFznlMzUp9YQcHDFUy7+wbu+s9lPl3LF3jSdn0JT/n1htAASHppJdm/XknOnzZwuERLV7HtCpyZn82SvhUErTZwQsNaU/agyn2w+vak1jh8gViQh2HXXjVwra09oKHbZgOcZlyhmg2sp+NZCPaFowC+/0+o+e9IPpr4d+Y/eYea+c0Ke56cSv6Nq+0PClDIvzBpNMqG7ZbnQpuWiv6wfzV1vMVJp3GYXXo9gaP19a7Zjtf1FZ0WMWkQRy8dB4oRcdQIt+UNzRrNsStLObq0DHFwkg0/czxMx6ljEMLDacqPR4yLo21hsdWFgyc0vGk7NMF/RLzi14sQC/IIW57ANXsvZOhb3RETZn6aDoE9f5tqbz8L0GBjiAqzeS48ERrBart+q3F4s6picee2GtXFuDiUzk7LEpw2KxP9wSrqPhhB444kcu9Ql7ScJpwxjSpiXBzHzxxF/Efb6Jg6Am2LHkFRENZsVVdVtFrEwUnoa45gmDkRfaRI+MemFRRBQIyNdWis6004s7z3RtJgXxAqqz7mdARHf1lGytPd2oagERAiIlm47gBfHyvkxPRjvSL4fUkzObCq4ganXVXificTHCUHNjQ22rwk+urDACSff4DcO1d1b7cSGuKQFIv20HX6BAynTaT85jFo9Ap0daH7agPCys2wequFn6LXg051i9a0G4g60Kgm/jEJsYPXjvHyzv2HXduZRsxQEBrgXd/2JszOg9ZCA1R+xtZW3l84lYa7slwKDUGrdZvuzxl8STMZrLbrtxrHeUOuVBprm93u523wmKuUdkKxmtdCWbsVNCK7Xyxi4fhNbLx3IlFVzQjNbRj2lgMQlxyDM34d84sJ/0TVOtrPnkLkpxuCmsDYFbdQwMnEr2N+MZouIxGHmzFs2xVwLprYWARBsGitgW67k07jKJjueKXEkrPCNIoKone3qK85ghgX5/A3Ze1WNOWHEBMTqfztVMJjOvjyrSlob67hwFkJCHoDYnIyglbrmJ+Jk1loAER8uAYhMpIDf3BTQ9YDi7pbCIJzbiGEk4lfxNdbCP9hG5ULHGaK8BtKZydEdqe8DFbb9VvBUb3dsaFIHJyk1l41V11zo34LujBLoJKgC3M5r2w7ZwrGrFSOvJrMUukjoiI6yZ5XTlpUI21Du8Bg5ODl+YiZ6Tb8xCEpaFOHIISpRlpz4B2AZkwhe58fTkeaHm32UMQC5+UfbRL59LhndxCKx4KiUHnnFA5VqEuFnuRT7Qs469tQgTf8lI4OjO3tZDzQPd2xDowD3NcKdpFwSOnosInyDVbb9VvBkZDmWCvQ1xxx7X3XY+RW9F2Wqu2dM8ai6PU2+Tk14wotn7Pv3MneS2JJuVvDjtY0rs5bwbZdmaxcVchFU1ZTcelQRp6zC+XYcWLmqWHyglaLsf4YR+fnWrxI22aOBqD13KnUTUkkLekEiRu0NExNZ/DLtqHeFj8QFyUbWqe4H2U0+6tBI5LzygHiE02rTCFa38VZ34YK/OVnft7MaP1XlJM9VQha16kDrBGstnM8hPUDdHU4f+i9SuNvZeMJ+/4nzN80UVG0zRiNxqCg2wLGaRNY92EUf1ryH+5JOZcDyyeyZ2YymZ9qiKpu4b/108heXseepgKafquQVg/7HiolruAYKX8MY9A/1VUaTXQ0UXuOYQCOXdrMsOuO0ro4mtgqA9E7atn01hgyh+zFcLQWFMUj24dllcYFDPXHELRa9FXVGPJVYWiupxtqcNW3oYBA8jNOm4DujwJqzm/H8MZGF6y267caR0ezj9myXKyxWwcXGVtbiTzcQuT2w5y4tIS9l+tYd/3j5IcdYePpTxNXdpTymsEcvqADYfNuhn7ehGH7bsQOhfx/NaDdWo8h2sCJ/YkcL1D9QARdGJqEeAyD1O/Drj0MHR2E3RdPc7pIS2Ey6d+csE0w42elepv7MwmhtvrQTUQDfvRtkBBIfu3J4Wi+3xiw81lzq7mlLGDntYOiKP3y74ILLri2rzn0V36hzG2AX//g1m81DtRSC6GMUOYXytxggJ8/CAq3/iw4BjCAAfQRBgTHAAYwAK/RnwXHC31NwA1CmV8oc4MBfv4gKNz6rcv5AAYwgL5Df9Y4BjCAAfQRBgTHAAYwAK/RLz1HJUmaBzwBiMCLsiw/EIRrZqEWy04FjMALsiw/IUnSIOBNIBvV/U+SZbnBdMzdwNWAAbhJluXPTNsnAS8DkaglMpd5UoTKA44isA6olmX57BDjloBaV2cMoABXAbtCgZ8kSbcA15h4bQWuBKL6ipskSf8AzgaOWpVdDVhfSpIUjvosTwLqgQtlWa7whmO/0zhML8fTwHxgFHCRJEmjgnBpPXCbLMsjgRLgl6br3gV8JctyPvCV6Tum3xYDo4F5wDMm7gDPoq6355v+5gWI4zJgh9X3UOL2BPCpLMuFwHgTzz7nJ0lSBnATMNn0koqma/clt5cdHBtIPlcDDbIs5wGPAX/1lmC/ExzAFGCvLMv7ZVnuBN4AFvb2RWVZPizL8gbT5ybUBz/DdO1XTLu9Apxj+rwQeEOW5Q5ZlstRK9hNMdXajZNleaVpNHrV6hifIUlSJnAW6qhuRqhwiwNORa3ahyzLnbIsHw8Vfqiad6SpdnEUal3jPuMmy/L/sK+fHEg+1ud6G5hlKs3qMfqj4MgADlp9rzJtCxokScoGioDVwBBTbVxM/80x0s54Zpg+99zuLx4H7kCdRpkRKtxygVrgn5IkbZQk6UVJkqJDgZ8sy9XAw0AlcBg4Icvy56HArQcCycdyjCzLeuAE4FXCkP4oOBxJxqCtKUuSFAO8A9wsy7KraDFnPAPOX5Ik83x4vYeHBI2bCVpgIvCsLMtFQAuuawQHs+0SUUfgHCAdiJYk6dJQ4OYhfOHjN9f+KDiqgCyr75moqmWvQ5IkHarQeF2W5XdNm4+Y1EJM/82hrc54Vpk+99zuD04BfiZJUgXq1O10SZJeCxFu5utVybJsriXwNqogCQV+s4FyWZZrZVnuAt4FykKEmzUCycdyjGl6Fo/91Mgl+qPgWAvkS5KUI0lSGKph6P3evqhpDvgSsEOW5UetfnofuNz0+XJgudX2xZIkhUuSlINqnFpjUjObJEkqMZ1zidUxPkGW5btlWc6UZTkbtT2+lmX50lDgZuJXAxyUJMlc2WgWsD1E+FUCJZIkRZnOOQvVfhUK3KwRSD7W5zof9Xk5uTUO05zsBuAz1A6WZVneFoRLnwJchjqabzL9nQk8AMyRJGkPMMf0HRMnGfUF+RT4pSzL5mQg16MaMfcC+4BPeolzKHG7EXhdkqQtwATgL6HAz6QFvQ1sQF2K1aC6bfcZN0mS/gOsBAokSaqSJOnqAPN5CUiSJGkvcCuup40OMeByPoABDMBr9DuNYwADGEDfY0BwDGAAA/AaA4JjAAMYgNcYEBwDGMAAvMaA4BjAAAbgNQYExwAGMACvMSA4BjCAAXiN/wMiJoNDW5TBBAAAAABJRU5ErkJggg==\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"from sentinel_helpers import scihub_cloud_mask\n",
|
|
"import matplotlib.pyplot as plt\n",
|
|
"\n",
|
|
"# pixels with clouds are True, pixels without are False\n",
|
|
"raster_cloud_mask, _ = scihub_cloud_mask(product_path)\n",
|
|
"plt.imshow(raster_cloud_mask[0])"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "markdown",
|
|
"metadata": {},
|
|
"source": [
|
|
"## Index Calculation"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 17,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"import numpy as np\n",
|
|
"import numpy.ma as ma\n",
|
|
"\n",
|
|
"import rasterio as r\n",
|
|
"from rasterio.enums import Resampling\n",
|
|
"import rasterio.plot as rplt\n",
|
|
"\n",
|
|
"from sentinel_helpers import RasterReaderList, scihub_cloud_mask\n",
|
|
"from tqdm.notebook import tqdm\n",
|
|
"\n",
|
|
"from tempfile import NamedTemporaryFile\n",
|
|
"\n",
|
|
"out_dir = Path('resources/spectral_indices/')\n",
|
|
"out_dir.mkdir(exist_ok=True, parents=True)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 18,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"application/vnd.jupyter.widget-view+json": {
|
|
"model_id": "ef37732b4a184e959bf49763e1e583aa",
|
|
"version_major": 2,
|
|
"version_minor": 0
|
|
},
|
|
"text/plain": [
|
|
"HBox(children=(HTML(value='Resampling datasets'), FloatProgress(value=0.0, max=2.0), HTML(value='')))"
|
|
]
|
|
},
|
|
"metadata": {},
|
|
"output_type": "display_data"
|
|
},
|
|
{
|
|
"name": "stdout",
|
|
"output_type": "stream",
|
|
"text": [
|
|
"\n",
|
|
"Wrote result to resources/spectral_indices/T33UUT_20180906T101021_NBR_10m.tif\n",
|
|
"CPU times: user 51.3 s, sys: 18.8 s, total: 1min 10s\n",
|
|
"Wall time: 39.6 s\n"
|
|
]
|
|
}
|
|
],
|
|
"source": [
|
|
"%%time\n",
|
|
"\n",
|
|
"with RasterReaderList(highest_resolution_band_paths) as readers:\n",
|
|
" # build the band_map as described above and scale up where needed\n",
|
|
" band_map = {}\n",
|
|
" for reader in tqdm(readers, desc='Resampling datasets'):\n",
|
|
" band_num = reader.name.split('_')[-2]\n",
|
|
" data_resolution = resolution(reader.name)\n",
|
|
" scale_factor = int(data_resolution / target_resolution)\n",
|
|
" out_shape = (\n",
|
|
" int(reader.height * scale_factor),\n",
|
|
" int(reader.width * scale_factor)\n",
|
|
" )\n",
|
|
" \n",
|
|
" band_map[band_num] = (np.clip(\n",
|
|
" # we read only the first band to obtain a two-dimensional ndarray\n",
|
|
" reader.read(1, out_shape=out_shape, resampling=Resampling.bilinear, masked=True),\n",
|
|
" 0, 10000) / 10000).astype('float32')\n",
|
|
" \n",
|
|
" if scale_factor == 1:\n",
|
|
" # this is the target resolution which defines the shape of the\n",
|
|
" # raster file the data is written to\n",
|
|
" out_name = Path(reader.name).name.replace(band_num, index_to_calculate.upper())\n",
|
|
" out_name = out_name.replace('.jp2', '.tif')\n",
|
|
" out_path = out_dir / out_name\n",
|
|
" out_meta = reader.meta.copy()\n",
|
|
" \n",
|
|
" # ignore numpy division errors (i.e. divide by 0) in this context;\n",
|
|
" # divide by zero results in `nan`. clouds are masked using the\n",
|
|
" # raster_cloud_mask created above\n",
|
|
" with np.errstate(divide='ignore', invalid='ignore'):\n",
|
|
" index = evaluate_formula(band_map, indices[index_to_calculate])\n",
|
|
" # we need to invert the cloud mask (`~`) because we want to hide those pixels that are cloudy\n",
|
|
" tile_mask = ~raster_cloud_mask | index.mask\n",
|
|
" index = ma.masked_array(index, tile_mask)\n",
|
|
" \n",
|
|
" out_meta.update({\n",
|
|
" 'count': 1,\n",
|
|
" 'driver': 'COG',\n",
|
|
" 'dtype': 'float32'\n",
|
|
" })\n",
|
|
" \n",
|
|
" with r.open(out_path, 'w+', **out_meta) as dst:\n",
|
|
" dst.write(index.data, 1)\n",
|
|
" # note that the mask we write has an inverted interpretation of \n",
|
|
" # True / False due to the `.msk` file format\n",
|
|
" dst.write_mask(~index.mask)\n",
|
|
" print(f'Wrote result to {out_path}')"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 19,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# to debug the cloud / tile mask, uncomment these lines\n",
|
|
"# plt.imshow(tile_mask)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 20,
|
|
"metadata": {},
|
|
"outputs": [
|
|
{
|
|
"data": {
|
|
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAQUAAAEDCAYAAAA8zxGMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfW0lEQVR4nO3df7xcdX3n8VckooDItRUwgvsgFjQgCkYWY33IaUQxrCCN5Xz4FVbSYhYoaiNGJZWfUXQ3stS2JMBGk5YY8IMSfljAsJUcHvZhCMqShWoUlrBKg6VsuLsRs9Lg7B/nzGbudebMueHe+b4Dn9fjMY9778yZOc+Zm/vJmZ9nUqvVIoqiqN3LUgOiKNIqhkIURSOKoRBF0YhiKERRNKIYClEUjSiGQhRFI5qcGtArM/sacALwlLsf3mB5Ay4FWsAGdz99YoVR9OJMeUthBTCryYJmdghwIfBud38L8GcTx4qiF3eyWwrufq+ZHdR5nJn9HnA1sC/wK+Cj7r4R+Chwtbs/U533qQFzo+hFk/KWQreuAz7m7u8APgUsqY5/E/AmM/sHM1tnZo22MKIo+u12maFgZq8Cfh+4ycweBK4FplQnTwYOAf4AOA1YZmZDg1dG0a6f7N2HLr0MGHb3I7uc9gSwzt3/FdhkZj+hHBL3D9AXRS+KdpktBXf/P5R/8DmAmU0ysyOqk28BZlbHv5by7sRjKZxRtKs3qcm7JM3scWAr8Dyw3d2PGnX6PsBK4N9Qbn182d2XV6fNAr4C7AYsc/cvNYGZ2Q2UdwdeC/wzcAnwXWAp5d2GlwM3uvvlZjYJuJLy2YrngS+4+41N1hNF0cjGMhSOcvene5y+ENjH3T9jZvsCPwFeR/kH+lPg/ZSb+PcDp7n7j8aHH0XReDdedx9awN7V/9ivArYA24GjgUfd/TF3fw64EThpnNYZRdEE1PSBxhawxsxawLXuft2o0/8auA3YDOwNnOLuvzGzA4Cfdyz3BPDObisws3nAPIBp06b9+Nhjj/0vwBCwB7AJmAGsBs4DLgKWZVl2VlEU12dZdmZRFCuAs4FFlE9VzgbWAVOBbcAwMA24Bzid8u7G4izLzi2KYmWWZXM6vl4FXAHMBe4EpgNPVtQpwAPA8cByYGGWZfO7XMZSYAFwAbCK8jGPjXGd4jopXKcsy9bTo6Z3H17v7pvNbD/gbsrXCtzbcfrJwLuBTwK/Vy1zBPAB4APufna13JnA0e7+sbr1FUXx8SzL/rIvLEHKNtD2KdsgfO0a3X1w983V16coJ+bRoxaZC9zs7i13f5Ryuk6j3DJ4Q8dyB1JuTfRrXRNXopRtoO1TtkH4gAZDwcz2MrO9298DxwEPj1rsZ8Cx1TL7A2+mfErwfuAQM5tqZrsDp1LezejX1MbXYPAp20Dbp2yD8AHNHlPYH1hdvgmRycAqd7/LzM4BcPdrKO9LrTCzh4BJwGfaz1SY2fnAdyifkvyau/9jg3VuG/M1GVzKNtD2KdsgfECDoeDuj1E+PjD6+Gs6vt9MuQXR7fx3AHeM0TU8xuUH2XBqQJ+GUwNqGk4N6NNwakCfhgexEtVXNE5LDahJ2QbaPmUbhA/QHQr3pAbUpGwDbZ+yDcIH6A4F5U9NUraBtk/ZBuEDdIfClakBNSnbQNunbIPwAbpDYXFqQE3KNtD2KdsgfEDDVzRGUfTSSXJLoSiKlakNvVK2gbZP2QbhaxdbClEUjSi2FMaYsg20fco2CF+72FKIomhEqlsKV6U29ErZBto+ZRuEr53kUKD8kAnVlG2g7VO2QfgA3aEwNzWgJmUbaPuUbRA+QHco3JkaUJOyDbR9yjYIH6A7FKanBtSkbANtn7INwgfoDoUn+y+SLGUbaPuUbRA+QHcoRFGUKNWhMKX/IslStoG2T9kG4QN0h8IDqQE1KdtA26dsg/ABukPh+NSAmpRtoO1TtkH4AN2hsDw1oCZlG2j7lG0QPkB3KCxMDahJ2QbaPmUbhA+IN0RFUTQqyS0F5bewKttA26dsg/C1iy2FKIpGFFsKY0zZBto+ZRuEr11sKURRNCLVLYWlqQ29UraBtk/ZBuFrJzkUgAWpATUp20Dbp2yD8AG6Q+GC1ICalG2g7VO2QfgA3aGwKjWgJmUbaPuUbRA+ACY3WcjMHge2As8D2939qFGnLwDO6LjMQ4F93X2Lmc0HzgZawEPAXHf/v31WORN4pOmVGHDKNtD2KdsgfEDDoVA1092f7naCuy+m2s+dmZ0IzK8GwgHAx4HD3H2bmTlwKrCiz7o2jsE16JRtoO1TtkH4gIm5+3AacEPHz5OBPcxsMrAnsLnBZQxNgGu8GkoN6NNQakBNQ6kBfRpKDejT0CBW0uh1Cma2CXiG8i7Ate5+XY/l9gSeAA529y3VcZ8AvgBsA9a4+xk9zjsPmAcwffr0X73rXe/6HOWNsAewCZgBrAbOAy4ClmVZdlZRFNdnWXZmURQrKO+mLAKWALOBdcDUat3DwDTgHuB0yt16L86y7NyiKFZmWTan4+tVlB+nPZfywzKns+OjsE4CrqF8G+tyYGGWZfO7XMZSykeLL6C8LziTctJP9HW6Ocuy48Z4naZQvld/Qq8T8O1qHYP4Pe3MdRquDIP4PSX9t5dl2Xp61Wq1+h7yPH999XW/PM835Hl+TI/lTsnz/PaOn1+T5/l38zzfN8/zl+d5fkue53P6rW/t2rVHN3GlOCjb1H3KtvDtODS6++Dum6uvT1UT8+gei57KyLsO7wM2ufu/uPu/AjcDv99glTOauBKlbANtn7INwgc0eKDRzPYCXubuW6vvjwMu77LcPkAGzOk4+mfAjOpuxTbgWOAHDVyrGyyTKmUbaPuUbRA+oNkDjfsD3zOzDcB64O/c/S4zO8fMzulYbjblYwbPto9w9/uAb1Let3uoWl/XxyNGdV7TK5AgZRto+5RtEL6y1PeTetx3mpzaMFbbe+7+VOuQmy5rqfoUDsq28O04qL6icVlqQE1dbfe+b/Gkn5588aTO4w72RSOe2jnslksH8ZbUXe62Eyp8xFunB9pbb7u49dCHLp/Uf8koSpfklkL1nLZkL8Q2eiC8+VuXj/tEfrHedhPdu9Z8ptX2vWvNZyT/pxzU7RdbComadvNlLYCNH74kthwiqVS3FFakNvRqvGwbP3zJpIkYCC+F226iCl+Z5FCgfMmoaso20PYp2yB8gO5QWJQaUJOyDbR9yjYIH6A7FJakBtSkbIPKd8hNixQfLNolbjvhBuJTHQqzUwNqUrZB5Xskv0jxAcxd4rYTbiA+1aGwLjWgJmUbjJPv6LsunIgtjZfEbTeBDcSnOhSmplz5wd/4/Ig/iLfceknnz0ltDRoX3/pZX5yILY2XxG03gQ3EpzoUtqVc+aOnfG4SlC8umnbzZa29XvHc/39dAYltDVL2KdsgfMDYPqNxkA2nBgD85I8u7va/5fCgHWNsODWgpuHUgD4Npwb0aXgQK1HdUpiWGlCTsg20fco2CB+gOxTuSQ2oSdkG2j5lG4QP0B0Kp6cG1KRsA22fsg3CB+gOhStTA2qqtR12y6WtqV+/IuULh3bZ206g8KE7FBanBtRUa/vRH146adMZC1O+cGiXve0ECh/x1ukoikYluaVQFMXK1IZeKdtA26dsg/C1iy2FAXTE7Re1Npy4aNLht13cejg+ji0SL7YUxtjO2DacuGjSEbdfNPqVkRPSi+22G2ThK4sthRfYm755eWv0pzj3KrYUol0h1S2Fq1IbejXa1nQgAAxiIOxKt51a4SuTHAqUe91VTdkG2j5lG4QP0B0Kc1MDalK2gbZP2QbhA3SHwp2pATUp20Dbp2yD8AG6Q2F6akBNyjbQ9inbIHyA7lB4MjWgJmUbaPuUbRA+QHcoRFGUqEafvGRmjwNbgeeB7e5+1KjTFwBndFzmocC+7r7FzIYo95Z7ONAC/tjdv99nlVOaXoEEKdtA26dsg/ABY/s4tpnu/nS3E9x9MdU7uMzsRGC+u2+pTv4KcJe7n2xmuwN7NljXA2NwDTplG2j7lG0QPmBiPqPxNOAGADN7NXAMcBaAuz8HPNfgMo4HHpoA23ikbANtn7INwgc0HwotYI2ZtYBr3f26bguZ2Z7ALOD86qg3Av8CLDezI4AfAp9w92f7rG95Q1eKlG2g7VO2QfiA5kPh3e6+2cz2A+42s43ufm+X5U4E/qHjrsNkyqdRPubu95nZV4DPAheNPqOZzQPmAcyYMWMv4D8AQ8AewCZgBrAaOK86/7Isy84qiuL6LMvOrPbIezbl/vaWUO5NZx3lZ+Vvo/wk3GmUn3N3OuWn2CzOsuzcoihWZlk2p+PrVZSvHptL+dzwdHY88vvJoigWUE7t5cDCLMvmd7mMpcAC4AJgFTAT2DiA6/RfgSPGeJ2mUG6aTuh1Ap4uiuI1A/o97cx1Oqgoiu8P6PeU9N9elmXr6dGY3xBlZpcCv3T3L3c5bTVwk7uvqn5+HbDO3Q+qfn4P8Fl3/+CYVhpF0cDq+5Skme1lZnu3vweOAx7ustw+QAbc2j7O3X8B/NzM3lwddSzwo37rVH4Lq7INtH3KNghfu75bCmb2RspNJyjvDqxy9y+Y2TkA7n5NtdxZwCx3P3XU+Y+kfEpyd+AxYK67PzOO1yGKovGs1WrJHdauXbsytWFXtKn7lG3h23GID1mJomhEki9zrh49lUzZBto+ZRuEr53kUKB8OkU1ZRto+5RtED5AdyhckBpQk7INtH3KNggfoDsUVqUG1KRsA22fsg3CB+gOhZmpATUp20Dbp2yD8AG6Q2FjakBNyjbQ9inbIHyA7lAYSg2oaSg1oE9D432Bb/7W5eP1vPXQOF3ORDWUGtCnoUGsRHUo7JEaUJOyDSbA95M/ar5viz695G67cW4gPtWhsCk1oCZlG2j7lG0QPkB3KMxIDahJ2QbaPmUbhA/QHQqr+y+SLGUbaPuUbRA+QHconJcaUJOyDbR9yjYIH6A7FH7rk5mEUraBtk/ZBuEDdIfCstSAmpRtoO1TtkH4gJ34OLYoil7cSW4pVB/yKZmyDbR9yjYIX7vYUoiiaESqWworUht6pWwDbZ+yDcLXTnIoUH6GvmrKNtD2KdsgfIDuUFiUGlCTsg20fco2CB+gOxSWpAbUpGwDbZ+yDcIH6A6F2akBNSnbQNunbIPwAbpDYV1qQE3KNtD2KdsgfIDuUJiaGlCTsg20fco2CB+gOxS2pQbUpGwDbZ+yDcIH6A6F4dSAmoZTA/o0nBpQ03BqQJ+GUwP6NDyIlagOhWmpATUp20Dbp2yD8AG6Q+Ge1ICalG2g7VO2QfgA3aFwempATco20PYp2yB8gO5QuDI1oCZlG2j7lG0QPgAmN1nIzB4HtgLPA9vd/ahRpy8Azui4zEOBfd19S3X6bsAPgH9y9xMarHIxcG4TW4KUbaDtU7ZB+ICGb52uhsJR7v50g2VPBOa7+3s7jvskcBTw6oZDIYqiRE3E3YfTgBvaP5jZgcAHGcNHSRVFsXICXOOSsg20fco2CF+7RncfgBawxsxawLXufl23hcxsT2AWcH7H0X8BfBrYu24FZjYPmAcwbdq0HwPHUO4maw/KnWDMoPyI6/MoP8ByWZZlZxVFcX2WZWdW7zU/m/KdZEsoXye+jvJVYNson+OdRvkI7umU988WZ1l2blEUK7Msm9Px9SrgCmAucCcwHXiyot5dFMVbgeOB5cDCLMvmd7mMpcACyt2Hr6LcOejGAVynrVD+AxrDdZoCPDCA6zRngL+nnblOlxVFMW9Av6ek//ayLFtPr1qtVt9Dnuevr77ul+f5hjzPj+mx3Cl5nt/e8fMJeZ4vqb7/gzzPv91kfWvXrl3ZZLlWq8Vht1zcarrseBzGYktxUPYp28K34zDmj2Mzs0uBX7r7l7ucthq4yd1XVT9/ETgT2A68Eng1cLO7zxnTSms62Be1HrWLxmtfh9E4deS3P9d6/jcv46EPXR6/m12svo8pmNleZrZ3+3vgOODhLsvtA2TAre3j3P1Cdz/Q3Q8CTgW+22QgVJtQjRr0QBiLLUUqvgdP+Pyk0QNBxdar8JU1eUxhf2C1mbWXX+Xud5nZOQDufk213Gxgjbs/Ow6uK8bhMiYqZRto+5RtEL6y1PeTetx3+nRqwyE3XdZStanfdruiLXw7DqqvaLwzNeCnJ1/c625JcluflH3KNggfoPsy5+mpATUp20Dbp2yD8AG6Q+HJ/ovU987vfLZ11J0LJ2JPNy/YNsEp+5RtED5Adyi84J799e784Pgr4umwAXbwNz4fuxt7EaQ6FKa80At4eOKeH3/Btgkume/RUz7X7zaP2+6FNRCf6lB4IDWgJmUbaPuUbRA+QHcoHJ8aUJOyDbR9yjYIH6A7FJanBtSkbANtn7INwgfoDoWFqQE1KdtA26dsg/ABDT9kJYqil06SWwrKH3ahbANtn7INwtcuthSiKBpRbCmMMWUbaPuUbRC+drGlEEXRiFS3FJamNvRK2QbaPmUbhK+d5FCg/NBJ1ZRtoO1TtkH4AN2hcEFqQE3KNtD2KdsgfIDuUFiVGlCTsg20fco2CB+gOxRmpgbUpGwDbZ+yDcIH6A6FjakBNSnbQNunbIPwAbpDYSg1oKah1IA+DaUG1DSUGtCnodSAPg0NYiWqQ2GP1ICalG2g7VO2QfgA3aGwKTWgJmUbaPuUbRA+QHcozEgNqEnZBto+ZRuED9AdCqtTA2pStoG2T9kG4QN0h8J5qQE1KdtA26dsg/ABukPhotSAmpRtoO1TtkH4AN2hsCw1oCZlG2j7lG0QPiDeOh1F0agktxSKorg+taFXyjbQ9inbIHztYkshiqIRTW6ykJk9DmwFnge2u/tRo05fAJzRcZmHAvsCewF/C7wO+A1wnbt/pd/6iqJYkWXZWc2uwmBTtoG2T9kG4Ws3lrsPM939yNEDAcDdF1enHQlcCBTuvgXYDlzg7odSvvDiT83ssAbrOnsMrkGnbANtn7INwgdMzGMKpwE3ALj7k+7+QPX9VuDHwAENLmPRBLjGK2UbaPuUbRA+oOHdB6AFrDGzFnCtu1/XbSEz2xOYBZzf5bSDgLcD9/U47zxgHsDhhx/+P4FjKN8Vtgfla75nUL6i6zzK52uXZVl2VlEU12dZdmZRFCsoJ+kiYAkwG1gHTAW2AcPANOAe4HTgSmBxlmXnFkWxMsuyOR1frwKuAOYCdwLTgScr6lNFUbyVcr9+y4GFWZbN73IZSyk/PusCyg/HmEn51teJvk4HQPnJv2O4TlMod146odcJWDjA39POXKc7iqKYN6DfU9J/e1mWradXrVar7yHP89dXX/fL83xDnufH9FjulDzPb+9y/KvyPP9hnucfbrK+tWvXfrzJcikOyjZ1n7ItfDsOje4+uPvm6utT1cQ8useip1LddWhnZi8HvgV83d1vbrI+yimrmrINtH3KNggf0OAxBTPby8z2bn8PHAc83GW5fYAMuLXjuEnAV4Efu/t/HoNr6hiWHXTKNtD2KdsgfECzxxT2B1abWXv5Ve5+l5mdA+Du11TLzQbWuPuzHed9N3Am8JCZPVgdt9Dd7+izzm0N/SlStoG2T9kG4QMaDAV3fww4osvx14z6eQWwYtRx3wMm7YRreCfOM6iGUwP6NJwaUNNwakCfhlMD+jQ8iJVIvsyZ8pFa1ZRtoO1TtkH4AN2hcE9qQE3KNtD2KdsgfIDuUDg9NaAmZRto+5RtED5AdyhcmRpQk7INtH3KNggfoDsUFqcG1KRsA22fsg3CB8Rbp6MoGpXklkJRFCtTG3qlbANtn7INwtcuthSiKBpRbCmMMWUbaPuUbRC+drGlEEXRiFS3FK5KbeiVsg20fco2CF87yaFA+SETqinbQNunbIPwAbpDYW5qQE3KNtD2KdsgfIDuULgzNaAmZRto+5RtED5AdyhMTw2oSdkG2j5lG4QP0B0KT/ZfJFnKNtD2KdsgfIDuUIiiKFGqQ2FKakBNyjbQ9inbIHyA7lB4IDWgJmUbaPuUbRA+QHcoHJ8aUJOyDbR9yjYIH6A7FJanBtSkbANtn7INwgfoDoWFqQE1KdtA26dsg/AB8YaoKIpGJbmloPwWVmUbaPuUbRC+drGlEEXRiGJLYYwp20Dbp2yD8LWLLYUoikakuqWwNLWhV8o20PYp2yB87SSHArAgNaAmZRto+5RtED5AdyhckBpQk7INtH3KNggf0GBX9ABm9jiwFXge2O7uR406fQFwRsdlHgrs6+5bzGwW8BVgN2CZu3+pwSpXNeMnSdkG2j5lG4QPaDgUqma6+9PdTnD3xVS7tDKzE4H51UDYDbgaeD/wBHC/md3m7j/qty7gkTHYBpmyDbR9yjYIHzAxdx9OA26ovj8aeNTdH3P354AbgZMaXMbGCXCNV8o20PYp2yB8QPMthRawxsxawLXufl23hcxsT2AWcH511AHAzzsWeQJ4Z4/zzgPmAbz97W/fClwCDAF7AJuAGcBq4DzgImBZlmVnFUVxfZZlZxZFsQI4G1gELAFmA+uAqcA2YBiYBtxDuUvvK4HFWZadWxTFyizL5nR8vYryk3PnUn4u3nR2fOrNvyuK4hnKd6wtBxZmWTa/y2UspXxg6ALKzb6ZlL/Uib5OnwXuHeN1mkL5ttwJvU7ATQP8Pe3MdfpFURTTBvR7SvpvL8uy9fSq1Wr1PeR5/vrq6355nm/I8/yYHsudkuf57R0/53meL+v4+cw8z/+q3/rWrl17ShNXioOyTd2nbAvfjkOjuw/uvrn6+lQ1MY/useip7LjrAOWWwRs6fj4Q2NxglZuauBKlbANtn7INwgc0eEzBzPYys73b3wPHAQ93WW4fIANu7Tj6fuAQM5tqZrtTDo3bGrhmNFgmVco20PYp2yB8QLMHGvcHvmdmG4D1wN+5+11mdo6ZndOx3Gxgjbs/2z7C3bdTPr7wHeDH5VH+jw3WubrxNRh8yjbQ9inbIHxlqe8n9bjv9MXUhl3Rpu5TtoVvx0H1FY0XpQbUpGwDbZ+yDcIH6L7MeVlqQE3KNtD2KdsgfEC8dTqKolFJbilUL3SRTNkG2j5lG4SvXWwpRFE0ItUthRWpDb1StoG2T9kG4WsnORQoX0eumrINtH3KNggfoDsUFqUG1KRsA22fsg3CB+gOhSWpATUp20Dbp2yD8AG6Q2F2akBNyjbQ9inbIHyA7lBYlxpQk7INtH3KNggfoDsUpqYG1KRsA22fsg3CB+gOhW2pATUp20Dbp2yD8AG6Q2E4NaCm4dSAPg2nBtQ0nBrQp+HUgD4ND2IlkkPh6quvnpba0CtlG2j7lG0QvnaSQ4HqA1xFU7aBtk/ZBuEDdIdCFEWJiqEQRdGIVIdC1/1KiKRsA22fsg3CB8Rbp6MoGpXqlkIURYmKoRBF0YjGstfp2szslcC9wCuqy/2mu19iZr8DfAM4CHgcMHd/pjrPhcCfUO7i/uPu/p3q+HcAKyj3e3cH8Al3b5nZK4C/Bd4B/C/gFHd/vDrPR4DPVZzPu/vfNPQtBk4EngP+BzDX3YcH6etl6zj9U5R79d63vedvhduuOu1jlPv22E65T5BPq/jM7EjgGuCVle88d18/aF+1zG7AD4B/cvcTVP4uujWeWwq/Bt7r7kcARwKzzGwG5Q5P/97dDwH+vvoZMzuMco9Rb6HcKe2S6oYDWEr5nOwh1WFWdfyfAM+4+8HAVcB/rC7rdyh3SPtOyl3aXWJmr2nouxs43N3fBvwUuDCBr5cNM3sD8H7gZ+2FVW47M5tJuRfxt7n7W4AvK/mA/wRc5u5HAhdXP6fwAXyCcodI7VT+Ln6rcRsK7t5y919WP768OrQo/9G0p9PfAH9YfX8ScKO7/9rdNwGPAkeb2RTg1e7+fXdvUU7AzvO0L+ubwLFmNgn4AHC3u2+ppu3d7LjBan3uvqbakxWU70I7cNC+mtsOyl/ypzt+lrntgHOBL7n7r6vlnhLztYBXV8fvw479mA7UZ2YHAh9k5Ee0S/xddGtcH1Mws93M7EHgqQpzH7C/uz8JUH3dr1q8227qD6gOT3Q5fsR5qj/k/w38bs1lNfF19seUu/8euK+bzcw+RLm5uWGUU+W2exPwHjO7z8wKM/u3Yr4/Axab2c8pt2IuTOT7C8rB/puO42T+LkY3rkPB3Z+vNtUOpJxuh9csPqnLca2a43f2PI18ZvbnlPc7v57C18X2NuDPKTd7R6dy200GXkO549MFgFf/Q6n4zgXmu/sbgPnAV1/AunbKZ2YnAE+5+w+7LNetgd92o5uQZx+qB+rWUm6q/HO16UP1tb2J2Ws39U+wYxO+8/gR5zGzyZSbhFtqLquJr/1gzAnAGdWmWTJfh+0kyvfPbzCzx6vzPGBmr0tlG+WbVZ335mrzfT3l/4SvFfJ9BLi5OukmyvvVI9Y1AN+7gQ9Vv8Mbgfea2UoE/y7ajdtQMLN9zWyo+n4P4H3ARspdz3+kWuwj7NhV/W3AqWb2CjObSvnAyfpqU2pr9UDWJODfjzpP+7JOBr5b/RF/BzjOzF5TPZByXHVcX5+ZzQI+A3zI3X/VcZaB+XrY/pu77+fuB7n7QZS/4Onu/guV2w64BXhvdfybgN2Bp4V8m4GsWuy9wCMd6xqIz90vdPcDq9/hqdX55iDyd9Gt8dxSmALcY2b/Hbif8n7dt4EvAe83s0coH0X/EoCXu6R34EfAXcCfuvvz1WWdS/mgzKOUTxO27+d/FfhdM3sU+CTVI7buvoXyk27vrw6XV8c18f01sDdwt5k9aGbXJPD1snVN6Lb7GvBGM3uY8n/Bj1RbDSq+jwJXmtkG4Aqqdxkm8HVL5e/it4qXOUdRNKJ4RWMURSOKoRBF0YhiKERRNKIYClEUjSiGQhRFI4qhEEXRiGIoRFE0ov8HSAOm0cPIMewAAAAASUVORK5CYII=\n",
|
|
"text/plain": [
|
|
"<Figure size 432x288 with 1 Axes>"
|
|
]
|
|
},
|
|
"metadata": {
|
|
"needs_background": "light"
|
|
},
|
|
"output_type": "display_data"
|
|
}
|
|
],
|
|
"source": [
|
|
"with r.open(out_path, 'r') as src:\n",
|
|
" rplt.show(src)"
|
|
]
|
|
},
|
|
{
|
|
"cell_type": "code",
|
|
"execution_count": 21,
|
|
"metadata": {},
|
|
"outputs": [],
|
|
"source": [
|
|
"# writing the mask explicitly works:\n",
|
|
"# dst.write(index.data, 1)\n",
|
|
"# dst.write_mask(~index.mask)"
|
|
]
|
|
}
|
|
],
|
|
"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
|
|
}
|