Monitoring and Mapping Wildfires Using Satellite Data
Contents
Monitoring and Mapping Wildfires Using Satellite Data#
Content creators: Brittany Engle, Diana Cadillo
Content reviewers: Yuhan Douglas Rao, Abigail Bodner, Jenna Pearson, Chi Zhang, Ohad Zivan
Content editors: Zane Mitrevica, Natalie Steinemann, Jenna Pearson, Chi Zhang, Ohad Zivan
Production editors: Wesley Banfield, Jenna Pearson, Chi Zhang, Ohad Zivan
Our 2023 Sponsors: NASA TOPS and Google DeepMind
# @title Project Background
from ipywidgets import widgets
from IPython.display import YouTubeVideo
from IPython.display import IFrame
from IPython.display import display
class PlayVideo(IFrame):
def __init__(self, id, source, page=1, width=400, height=300, **kwargs):
self.id = id
if source == "Bilibili":
src = f"https://player.bilibili.com/player.html?bvid={id}&page={page}"
elif source == "Osf":
src = f"https://mfr.ca-1.osf.io/render?url=https://osf.io/download/{id}/?direct%26mode=render"
super(PlayVideo, self).__init__(src, width, height, **kwargs)
def display_videos(video_ids, W=400, H=300, fs=1):
tab_contents = []
for i, video_id in enumerate(video_ids):
out = widgets.Output()
with out:
if video_ids[i][0] == "Youtube":
video = YouTubeVideo(
id=video_ids[i][1], width=W, height=H, fs=fs, rel=0
)
print(f"Video available at https://youtube.com/watch?v={video.id}")
else:
video = PlayVideo(
id=video_ids[i][1],
source=video_ids[i][0],
width=W,
height=H,
fs=fs,
autoplay=False,
)
if video_ids[i][0] == "Bilibili":
print(
f"Video available at https://www.bilibili.com/video/{video.id}"
)
elif video_ids[i][0] == "Osf":
print(f"Video available at https://osf.io/{video.id}")
display(video)
tab_contents.append(out)
return tab_contents
video_ids = [('Youtube', '2CDVn6O_q34'), ('Bilibili', 'BV1GP411k7XA')]
tab_contents = display_videos(video_ids, W=730, H=410)
tabs = widgets.Tab()
tabs.children = tab_contents
for i in range(len(tab_contents)):
tabs.set_title(i, video_ids[i][0])
display(tabs)
# @title Tutorial slides
# @markdown These are the slides for the videos in all tutorials today
from IPython.display import IFrame
link_id = "8w5da"
In this project, you will be working with wildfire remote sensing data from Sentinel 2 to extract Burn Perimeters using multiple Burn Indices and other relevant information related to wildfires. By integrating this data with information from other global databases, you will have the opportunity to explore the connections between wildfires and their impact on both the ecosystem and society.
Project Template#
Note: The dashed boxes are socio-economic questions.
Data Exploration Notebook#
Project Setup#
# installs
# !pip install gdal
# !pip install pandas
# !pip install geopandas
# imports
import random
import numpy
import matplotlib.pyplot as plt
import pandas as pd
import geopandas as gpd
import os
import matplotlib
/tmp/ipykernel_419586/2929139239.py:6: DeprecationWarning: Shapely 2.0 is installed, but because PyGEOS is also installed, GeoPandas still uses PyGEOS by default. However, starting with version 0.14, the default will switch to Shapely. To force to use Shapely 2.0 now, you can either uninstall PyGEOS or set the environment variable USE_PYGEOS=0. You can do this before starting the Python process, or in your code before importing geopandas:
import os
os.environ['USE_PYGEOS'] = '0'
import geopandas
In the next release, GeoPandas will switch to using Shapely by default, even if PyGEOS is installed. If you only have PyGEOS installed to get speed-ups, this switch should be smooth. However, if you are using PyGEOS directly (calling PyGEOS functions on geometries from GeoPandas), this will then stop working and you are encouraged to migrate from PyGEOS to Shapely 2.0 (https://shapely.readthedocs.io/en/latest/migration_pygeos.html).
import geopandas as gpd
Satellite Data#
Global Wildfire Information System (GWIS) is a global, joint initiative created by the GEO and the Copernicus Work Programs. The goal of this program is to bring together sources at a regional and international level to provide a global, comprehensive understanding of fire regimes and their effects.
The Globfire dataset uses the MODIS burned area product (MCD64A1) to determine the amount of burnt area per event (Artés et al., 2019). The MCD64A1 product that Globfire is built on top of combines imagery from the Terra and Aqua satellites with thermal anomalies to provide burnt area information (Website, MODIS C6 BA User Guide & User Guide).
Sentinel 2 Launched by the European Space Agency in June of 2015 (S2-A) and March of 2017 (S2-B), the Copernicus Sentinel-2 mission combines two polar-orbiting satellites to monitor variability in land surface conditions. Together, they have a global revisit time of roughly 5 days.
In the provided Project Template, we load the following datasets:
Global Wildfire Information Systems: Climate Action Large Wildfire Dataset#
The Climate Action Large Wildfire Dataset is a filtered version of the Globfire dataset (created by GWIS, which you can access here).
The resolution of this dataset is 500m. It has been pre-filtered to include only fires that are Class F or greater. Per the National Wildfire Coordinating Group, a Class F fire is defined as a fire that is 1,000 acres or greater, but less than 5,000 acres (note that a 2000m square region is roughly 1000 acres). For this dataset, all fires greater than 1,000 acres are included. Additional columns indicating area (acres), landcover number and landcover description, and country in which they are located within, were added. The landcover number and description were added using the Copernicus Global Land Cover Layers: CGLS-LC100 Collection 3 User Guide.
Table Information: ID: Globfire fire ID (unique to each fire)
IDate: Globfire Initial (start) date of the fire
FDate: Globfire Final (end) date of the fire
Continent: Location of the fire
Area_acres: size of fire in acres
Landcover: land cover value from ESA, if the landcover of the fire cover is greater than 51%, then it is labeled as that landcover
LC_descrip: land cover description from ESA
Country: country the fire is located in
# code to retrieve and load the data
url_climateaction = "~/shared/Data/Projects/Wildfires/ClimateAction_countries.shp"
Dataset = gpd.read_file(url_climateaction) # need to update to OSF and pooch.retrieve
---------------------------------------------------------------------------
CPLE_OpenFailedError Traceback (most recent call last)
File fiona/ogrext.pyx:136, in fiona.ogrext.gdal_open_vector()
File fiona/_err.pyx:291, in fiona._err.exc_wrap_pointer()
CPLE_OpenFailedError: /home/wesley/shared/Data/Projects/Wildfires/ClimateAction_countries.shp: No such file or directory
During handling of the above exception, another exception occurred:
DriverError Traceback (most recent call last)
Cell In[5], line 3
1 # code to retrieve and load the data
2 url_climateaction = "~/shared/Data/Projects/Wildfires/ClimateAction_countries.shp"
----> 3 Dataset = gpd.read_file(url_climateaction) # need to update to OSF and pooch.retrieve
File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/geopandas/io/file.py:281, in _read_file(filename, bbox, mask, rows, engine, **kwargs)
278 else:
279 path_or_bytes = filename
--> 281 return _read_file_fiona(
282 path_or_bytes, from_bytes, bbox=bbox, mask=mask, rows=rows, **kwargs
283 )
285 else:
286 raise ValueError(f"unknown engine '{engine}'")
File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/geopandas/io/file.py:322, in _read_file_fiona(path_or_bytes, from_bytes, bbox, mask, rows, where, **kwargs)
319 reader = fiona.open
321 with fiona_env():
--> 322 with reader(path_or_bytes, **kwargs) as features:
323 crs = features.crs_wkt
324 # attempt to get EPSG code
File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/fiona/env.py:457, in ensure_env_with_credentials.<locals>.wrapper(*args, **kwds)
454 session = DummySession()
456 with env_ctor(session=session):
--> 457 return f(*args, **kwds)
File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/fiona/__init__.py:292, in open(fp, mode, driver, schema, crs, encoding, layer, vfs, enabled_drivers, crs_wkt, allow_unsupported_drivers, **kwargs)
289 path = parse_path(fp)
291 if mode in ("a", "r"):
--> 292 colxn = Collection(
293 path,
294 mode,
295 driver=driver,
296 encoding=encoding,
297 layer=layer,
298 enabled_drivers=enabled_drivers,
299 allow_unsupported_drivers=allow_unsupported_drivers,
300 **kwargs
301 )
302 elif mode == "w":
303 colxn = Collection(
304 path,
305 mode,
(...)
314 **kwargs
315 )
File ~/miniconda3/envs/climatematch/lib/python3.10/site-packages/fiona/collection.py:243, in Collection.__init__(self, path, mode, driver, schema, crs, encoding, layer, vsi, archive, enabled_drivers, crs_wkt, ignore_fields, ignore_geometry, include_fields, wkt_version, allow_unsupported_drivers, **kwargs)
241 if self.mode == "r":
242 self.session = Session()
--> 243 self.session.start(self, **kwargs)
244 elif self.mode in ("a", "w"):
245 self.session = WritingSession()
File fiona/ogrext.pyx:588, in fiona.ogrext.Session.start()
File fiona/ogrext.pyx:143, in fiona.ogrext.gdal_open_vector()
DriverError: /home/wesley/shared/Data/Projects/Wildfires/ClimateAction_countries.shp: No such file or directory
We can now visualize the content of the dataset.
# code to print the shape, array names, etc of the dataset
Dataset.head()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[6], line 2
1 # code to print the shape, array names, etc of the dataset
----> 2 Dataset.head()
NameError: name 'Dataset' is not defined
# plot the dataset
Dataset.plot()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 2
1 # plot the dataset
----> 2 Dataset.plot()
NameError: name 'Dataset' is not defined
Great work! Now you are all set to address the questions you are interested in! Good luck!
Sentinel-2#
Each of the two satellites contain a Multi-Spectral Instrument (MSI) that houses 13 spectral bands, at various resolutions. The MSI uses a push-broom acquisition technique and measures in the Visible and Near Infrared (NIR) and the Short Wave Infrared (SWIR) domains. There are:
Four bands with 10m resolution: B2-Blue (490 nm), B3- Green (560 nm), B4- Red (665 nm) and Band 8- Near Infrared (NIR) (842 nm).
Six bands with 20m resolution: B5, B6, B7 and B8A (705 nm, 740 nm, 783 nm and 865 nm respectively) Vegetation Red Edge bands, along with B11 and B12 (1610 nm and 2190 nm) SWIR bands. These bands are mostly used for vegetation characterisation, vegetation stress assessment and ice/cloud/snow detection.
Three additional bands: B1 - Coastal Aerosol (443 nm), B9- Water Vapor (945 nm), and B10-SWIR-Cirrus (1375 nm) which are typically used for corrections of the atmosphere and clouds.
# view image
def showImage(Output):
plt.imshow(Output)
plt.show()
# data source-specific imports
# root folder location of where the imagery is currently saved
rootFolder = '/home/jovyan/shared/Data/Projects/Wildfires'
continet = "Asia"
# import pre images
# pre_fire_paths = glob.glob(rootFolder + continet + id +"/pre_fire_*.npy")
pre_fire_paths = [
fileName
for fileName in os.listdir(os.path.join(rootFolder, continet))
if (fileName.endswith(".npy") & fileName.startswith("pre_fire_"))
]
pre_fires_numpy = [
numpy.load(os.path.join(rootFolder, continet, x)).astype(int)
for x in pre_fire_paths
]
# import post images
post_fire_paths = [
fileName
for fileName in os.listdir(os.path.join(rootFolder, continet))
if (fileName.endswith(".npy") & fileName.startswith("post_fire_"))
]
post_fires_numpy = [
numpy.load(os.path.join(rootFolder, continet, x)).astype(int)
for x in post_fire_paths
]
# import Pre-SCL Mask for cloud coverage
scl_fire_paths = [
fileName
for fileName in os.listdir(os.path.join(rootFolder, continet))
if (fileName.endswith(".npy") & fileName.startswith("scl_pre_fire_"))
]
scl_fires_numpy = [
numpy.load(os.path.join(rootFolder, continet, x)) for x in scl_fire_paths
]
# import Post-SCL Mask for cloud coverage
scl_fires_post_paths = [
fileName
for fileName in os.listdir(os.path.join(rootFolder, continet))
if (fileName.endswith(".npy") & fileName.startswith("scl_post_fire_"))
]
scl_fires_post_numpy = [
numpy.load(os.path.join(rootFolder, continet, x)) for x in scl_fires_post_paths
]
---------------------------------------------------------------------------
FileNotFoundError Traceback (most recent call last)
Cell In[9], line 12
6 continet = "Asia"
8 # import pre images
9 # pre_fire_paths = glob.glob(rootFolder + continet + id +"/pre_fire_*.npy")
10 pre_fire_paths = [
11 fileName
---> 12 for fileName in os.listdir(os.path.join(rootFolder, continet))
13 if (fileName.endswith(".npy") & fileName.startswith("pre_fire_"))
14 ]
15 pre_fires_numpy = [
16 numpy.load(os.path.join(rootFolder, continet, x)).astype(int)
17 for x in pre_fire_paths
18 ]
20 # import post images
FileNotFoundError: [Errno 2] No such file or directory: '/home/jovyan/shared/Data/Projects/Wildfires/Asia'
# print list of pre_fires
print("\n".join(pre_fire_paths))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[10], line 2
1 # print list of pre_fires
----> 2 print("\n".join(pre_fire_paths))
NameError: name 'pre_fire_paths' is not defined
# print list of post_fire
print("\n".join(post_fire_paths))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[11], line 2
1 # print list of post_fire
----> 2 print("\n".join(post_fire_paths))
NameError: name 'post_fire_paths' is not defined
# print list of SCL
print("\n".join(scl_fire_paths))
print("\n".join(scl_fires_post_paths))
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[12], line 2
1 # print list of SCL
----> 2 print("\n".join(scl_fire_paths))
3 print("\n".join(scl_fires_post_paths))
NameError: name 'scl_fire_paths' is not defined
# view pre-fire satellite image that was taken right before the fire start date
showImage(
numpy.clip(pre_fires_numpy[2][:, :, [3, 2, 1]] / 10000 * 3.5, 0, 1)
) # RGB bands for Sentinel 2 are Bands: 4,3,2
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[13], line 3
1 # view pre-fire satellite image that was taken right before the fire start date
2 showImage(
----> 3 numpy.clip(pre_fires_numpy[2][:, :, [3, 2, 1]] / 10000 * 3.5, 0, 1)
4 ) # RGB bands for Sentinel 2 are Bands: 4,3,2
NameError: name 'pre_fires_numpy' is not defined
# view post-fire satellite image that was taken right before the fire start date
showImage(
numpy.clip(post_fires_numpy[0][:, :, [3, 2, 1]] / 10000 * 3.5, 0, 1)
) # RGB bands for Sentinel 2 are Bands: 4,3,2
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[14], line 3
1 # view post-fire satellite image that was taken right before the fire start date
2 showImage(
----> 3 numpy.clip(post_fires_numpy[0][:, :, [3, 2, 1]] / 10000 * 3.5, 0, 1)
4 ) # RGB bands for Sentinel 2 are Bands: 4,3,2
NameError: name 'post_fires_numpy' is not defined
Sentinel-2: Vegetation Health & Burnt Areas#
Sentinel-2 Imagery - Vegetation Health & Burnt Areas
Continents included:
Asia
Africa
Austrailia
Europe
North America
South America
Vegetation Health: Vegetation indices define and monitor the health of vegetation using the radiance values of the visible and near-infrared channels. The Normalized Difference Vegetation Index (NDVI) is used to measure “greenness” of vegetation. As one of the most widely used vegetation indexes, the NDVI takes advantage of how strongly Chlorophyll absorbs visible light and how well the leafs cellular structure reflects near-infrared. Its values range from +1.0 to -1.0, with areas of sand, rock, and snow typically having a value of 0.1 or less. Shrubs and spare vegetation are roughly 0.2 to 0.5, while higher NDVI values (0.6 to 0.9) indicate dense vegetation that is typically found in tropical forests.
The NDVI can also be used to create the Vegetation Condition Index (VCI). The VCI depends on the current NDVI along with the extreme NDVI values within the dataset (NDVImax and NDVImin). Specifically, $\(VCI = \frac{NDVI-NDVImin}{NDVImax-NDVImin}\times 100\%\)$
This number is then compared to the threshold to determine the drought severity. For this project, up to 3 months of pre-fire imagery will be used to determine current drought conditions.
Burnt Area Severity: Burn severity describes how the intensity of the fire affects the functioning of the ecosystem in which it occurred. The degree to which it alters the ecosystem is typically found using a burn index, which then (typically) classes the severity of the fire as: unburned, low severity, moderate severity, or high severity. One of the most common burn indices is the Normalized Burn Ratio (NBR). This index is designed to highlight burnt areas in fire zones. The formula is similar to NDVI, except that the formula combines the use of both near infrared (NIR) and shortwave infrared (SWIR) wavelengths. Healthy vegetation shows a very high reflectance in the NIR, and low reflectance in the SWIR portion of the spectrum - the opposite of what is seen in areas devastated by fire. Recently burnt areas demonstrate low reflectance in the NIR and high reflectance in the SWIR. The difference between the spectral responses of healthy vegetation and burnt areas reach their peak in the NIR and the SWIR regions of the spectrum. The difference between normalized burn ratios before and after a fire is called the dNBR, and is one of many useful indices.
Specifically, the dNBR isolates the burned areas from the unburned areas, and subtracts the pre-fire imagery from the post-fire imagery. This removes any unchanged, and thus unburned, pixels as their values result in near zero. The results of the dNBR are based on burn severity, and correspond to the gradient of burn severity for every pixel. The dNBR has an unscaled range of -2.0 to +2.0 with burned areas tending to show more positively.
SCL Mask Importing the SCL in order to mask out clouds form imagery
# compute SCL Mask to 0s and 1s, masking out clouds and bad pixels
def computeSCLMask(image):
rImage, cImage = image.shape
sclOutput = numpy.zeros((rImage, cImage))
for x in range(cImage):
for y in range(rImage):
sclOutput[y, x] = 1 if image[y, x] in [0, 1, 3, 8, 9, 11] else 0
return sclOutput
# create Pre-fire and post-fire SCL masks
print(scl_fires_numpy[5])
pre_SCL_Mask = computeSCLMask(scl_fires_numpy[5])
post_SCL_Mask = computeSCLMask(scl_fires_post_numpy[0])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[16], line 2
1 # create Pre-fire and post-fire SCL masks
----> 2 print(scl_fires_numpy[5])
3 pre_SCL_Mask = computeSCLMask(scl_fires_numpy[5])
4 post_SCL_Mask = computeSCLMask(scl_fires_post_numpy[0])
NameError: name 'scl_fires_numpy' is not defined
# view SCL imagery for image closests to the start fire date
# ------- Save SCL as colored image based on SCL classification
# No Data (0) = black
# Saturated / Defective (1) = red
# Dark Area Pixels (2) = chocolate
# Cloud Shadows (3) = brown
# Vegetation (4) = lime
# Bare Soils (5) = yellow
# Water (6) = blue
# Clouds low probability / Unclassified (7) = aqua
# clouds medium probability (8) = darkgrey
# Clouds high probability (9) light grey
# Cirrus (10) = deepskyblue
# Snow / Ice (11) = magenta
# colors: https://matplotlib.org/3.1.1/gallery/color/named_colors.html#sphx-glr-gallery-color-named-colors-py
def showSCL(image):
cmap = matplotlib.colors.ListedColormap(
[
"black",
"red",
"chocolate",
"brown",
"lime",
"yellow",
"blue",
"aqua",
"darkgrey",
"lightgrey",
"deepskyblue",
"magenta",
]
)
plt.imshow(image, cmap=cmap)
plt.show()
showSCL(scl_fires_numpy[5])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[17], line 40
36 plt.imshow(image, cmap=cmap)
37 plt.show()
---> 40 showSCL(scl_fires_numpy[5])
NameError: name 'scl_fires_numpy' is not defined
Dataset-specific Functions: NDVI#
# this code computes the NDVI for an image
"The NDVI can be computed on any image (pre or post)."
"Compute the NDVI on the pre-fire image"
# compute NDVI
def computeNDVI(image, mask):
r, c, ch = image.shape
ndviOutput = numpy.zeros((r, c))
for x in range(c):
for y in range(r):
if (image[y, x, 7] == 0 and image[y, x, 3] == 0) or mask[y, x] == 1:
ndviOutput[y, x] = numpy.nan
else:
ndviOutput[y, x] = (image[y, x, 7] - image[y, x, 3]) / (
image[y, x, 7] + image[y, x, 3]
)
return ndviOutput
# TA Code
computeNDVI_value = computeNDVI(pre_fires_numpy[2], pre_SCL_Mask)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[19], line 2
1 # TA Code
----> 2 computeNDVI_value = computeNDVI(pre_fires_numpy[2], pre_SCL_Mask)
NameError: name 'pre_fires_numpy' is not defined
# plot NDVI without remap
showImage(computeNDVI_value)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[20], line 2
1 # plot NDVI without remap
----> 2 showImage(computeNDVI_value)
NameError: name 'computeNDVI_value' is not defined
# use this code to remap the NDVI to specific colors for values
def remapNDVI(NDVI):
remapped = numpy.zeros((NDVI.shape[0], NDVI.shape[1]))
for x in range(remapped.shape[0]):
for y in range(remapped.shape[1]):
if numpy.isnan(NDVI[x, y]):
remapped[x, y] = numpy.nan
elif NDVI[x, y] <= -0.2:
remapped[x, y] = 1
elif NDVI[x, y] <= 0:
remapped[x, y] = 2
elif NDVI[x, y] <= 0.1:
remapped[x, y] = 3
elif NDVI[x, y] <= 0.2:
remapped[x, y] = 4
elif NDVI[x, y] <= 0.3:
remapped[x, y] = 5
elif NDVI[x, y] <= 0.4:
remapped[x, y] = 6
elif NDVI[x, y] <= 0.5:
remapped[x, y] = 7
elif NDVI[x, y] <= 0.6:
remapped[x, y] = 8
elif NDVI[x, y] <= 0.7:
remapped[x, y] = 9
elif NDVI[x, y] <= 0.8:
remapped[x, y] = 10
elif NDVI[x, y] <= 0.9:
remapped[x, y] = 11
elif NDVI[x, y] <= 1:
remapped[x, y] = 12
else:
remapped[x, y] = 13
return remapped
# TA Code
NDVI_remap = remapNDVI(computeNDVI_value)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[22], line 2
1 # TA Code
----> 2 NDVI_remap = remapNDVI(computeNDVI_value)
NameError: name 'computeNDVI_value' is not defined
# view remapped NDVI
def showNDVI(image):
cmap = matplotlib.colors.ListedColormap(
[
"#000000",
"#a50026",
"#d73027",
"#f46d43",
"#fdae61",
"#fee08b",
"#ffffbf",
"#d9ef8b",
"#a6d96a",
"#66bd63",
"#1a9850",
"#006837",
]
)
plt.imshow(image, cmap=cmap)
plt.show()
showNDVI(NDVI_remap)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[23], line 23
19 plt.imshow(image, cmap=cmap)
20 plt.show()
---> 23 showNDVI(NDVI_remap)
NameError: name 'NDVI_remap' is not defined
Dataset-specific Functions: VCI#
# compute the SCL mask for all the SCLs and apply it to the pre_fire_NDVIs
pre_fires_scl = [computeSCLMask(x) for x in scl_fires_numpy]
pre_fires_NDVI = [computeNDVI(x[0], x[1]) for x in zip(pre_fires_numpy, pre_fires_scl)]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[24], line 2
1 # compute the SCL mask for all the SCLs and apply it to the pre_fire_NDVIs
----> 2 pre_fires_scl = [computeSCLMask(x) for x in scl_fires_numpy]
3 pre_fires_NDVI = [computeNDVI(x[0], x[1]) for x in zip(pre_fires_numpy, pre_fires_scl)]
NameError: name 'scl_fires_numpy' is not defined
# compute for VCI
def computeVCI(for_ndvi_image, ndvi_dataset):
rImage, cImage = for_ndvi_image.shape
vciOutput = numpy.zeros((rImage, cImage))
ndvi_dataset.append(for_ndvi_image)
for x in range(cImage):
for y in range(rImage):
pixels = [z[y, x] for z in ndvi_dataset]
filtered_pixels = [f for f in pixels if not numpy.isnan(f)]
if len(filtered_pixels) == 0 or len(filtered_pixels) == 1:
vciOutput[y, x] = 1
elif numpy.isnan(for_ndvi_image[y, x]):
vciOutput[y, x] = 1
else:
max_val = max(filtered_pixels)
min_val = min(filtered_pixels)
if max_val == min_val:
vciOutput[y, x] = 1
else:
vciOutput[y, x] = (for_ndvi_image[y, x] - min_val) / (
max_val - min_val
)
return vciOutput
VCI Drought Threshold#
image and more information: Spatial and Temporal Variation of Drought Based onSatellite Derived Vegetation Condition Index in Nepal from 1982–2015
![droughtthreshold.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAjcAAAD5CAYAAADbVRRDAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAALEoAACxKAXd6dE0AAEvrSURBVHhe7d1/aBt3nj/+Z770DhuynHJkuTG0EIUUViaByOyCZS5/eEKWs0IKkXEh8id7JEoLqdxCKyeQ2u1BKne5VM5C1m6XrpVwLVKgQSq0WOEarMBlkQJZpECCVWiwAglIR3pIkIAFDby/M5q3bEmWbPlXbCvPB7yTmdFoRrLmx2tm3u/Xe5vQgIiIiKhJ/H/yfyIiIqKmwOCGiIiImgqDGyIiImoqDG6IiIioqTC4ISIioqbC4IaIiIiaCoMbIiIiaioMboiIiKipMLghIiKipsLghoiIiJoKgxsiIiJqKgxuiIiIqKmse8eZT58+xf/+7//KMSIiImp2//zP/1wsG2Xdg5v//u//xn/+53/KMSIiImp2//7v/14sG2Xdg5v/+Z//wV/+8hc5RkRERM2ut7e3WDbKugc3z58/R6FQkGNERETU7P7xH/+xWDbKugc3RERERC8SW0sRERFRU2FwQ0RERE2FwQ0RERE1FQY3RERE1FQY3BAREVFTYXBDRERETYXBDRERETUVBjdERETUVBjcEBERUVNhcENERERNhcENERERNRUGNy+lAvI/RhG8mZbjRETNKf+3IKIP5cgayN+PIno3K8fqSSN6NY68HKMXjx1n1lF4mEQsmUAyqf1/O1W1kZphPWSBZY8VHdYOWHaZ0CJf2ZR+TiF6L4O89p0Sd6IIfxtBSt83z8cgPrIZ89CyFB5r28VPOTm2Ojte74L11U29BdFmJ/fxxe1Au80KZcGmlkfqZgI13/2rdnT9Vln8+FbQ3n83itgt7Xh5J45UxcHSBEuntn1bu9B1wArLzrIlNfSZ29DRbdGWshIFpK4M4PDJHIbSIbh2yclS/u9+eC/4Eb4W10IRBZYjvXC/Pwx3tyLnqOHnME792gtzPIahzkX+Ks+TGD3QgXHLBCYvuWDZLqfTi6MHN1SSE9Nhr3DuVfSAb74oFqEeUo1S/Zp83XkmIBJP5WI2mztjtT/7+ZicgZbt1lDl33IVZeiWXCbRSt2bEA5tH7ftqb2NmTv149eQmMzJ+StMi4k+/XWbMJe9R9mrTTs3qR0Va8vdCwnvMYvQQoGydSnC0q0vSyud5rLpRlH2OoX3hxljAfK4ZFEq55mfV1/OmEgYcy/TrEhctBc/m+Mrub4yM0Gn8bl7fMZx+5cZETiuHx8VoX4a095dW+KCVaAvIDJyfFE/TQi7/j2OyHXQC8XgRppNBYSnu+zkr6jCdWlSTD+psZn/MisyqSkxccZedjBQxcRP8vVNKyaG5j6vVhjcrNwdb/HgaDvhE4HIlJiKlpVvhoRa/nfWto2hb6rmiQSE97h2oNRe98blMolW6+mU8CwIFrzant+Y6c9txnt6JsTCkEDSAoHQ+2pFUGPu84nJVK5GUKAdK+MTlcfWM1NV802LsUPzy9KL9cLKQpqSueBlvxZYyGlzkj5hk+txf1cWusWNfVoPcIaiNY77T0LCBau2v9YLfarNith5Yx9Xjgfq/z1pXTC40cxqG3VpYy+WziEx9US+uJR0SLg79fcNNXwA2Tgx4S3/ngxuVu5RQDjqBbTF18r+znCIwCP5Wjntyk6t9xrRCs1ctpdte3qxN3jhNSumzujz1zm5654mhO9I+R1gm3CHGzlt5xa5k5IRgb7S8oziCDZ0b6S2suDFfnnhZ0t8Wvr8Vftv+X57uvqOlQxUGr1rUyLv3hQ/y5cMb14kVih+GITr6DDichSKE4GrXqg75fhSdjnQr229OGSGtuPSS8WE1tVUlWlpXWFdAqL6zMc8GKo4GEUQvtVA44F8FOHPtP/3D6K3u9aGnUbwHTsGvy9VplXgDAYwdtQsxxdjgrWvH13akFlpMyatiyyCnw7K47l2bNaOy5WySCVLn79q/33VjA45iHhK+7ZlHofh+xjwfehc3nF+z2H09xmDkbcHEXxsDNP6e8mDmzT87/QjOFfxXYH7y3E4qyqeLUVRVG0jbsMOOU4vgZ0qhqPDjQfBtazFMoiqtajo/cAqRwyRS2Ek5XA9+ZthjGv/299zoPLdhvz3o/B8XdZKqG8U3mONBDaSFjxYtCW3/Xo1VwRLuBuA75oc3t8Fy6tyuEGz8n/82lR24VFA/IoP4RPDcO2XkxqmwGIr/TW1AOnqUr8CrZWXOrgp3PRj+Loc0SkD6D+y/Gvpth4fpt7rWtlVeCGPwnM5XEchv/Q8dWnLzz+Tw6uxVstpFi0KrN21Wp4sw1osoxHPCyv/7Vb4u+vbLLeXjWM96oZ+Q3nOXR9CNwtypJY0Qn/1a/+74DpSK2BJwv/xOMpCG3jecmAZoY3Ggv6oH/2/kaNrroDoN775IM5mqfH5yoONNDI/y0Hd4zRSclD5rRlz95dKd220oG8lx3jzHv1+lSF5MYToYj8DrZmXOrhJ/s1fsbPiuK3mFctSWl61Qv1N/c2+kE8j+a0fgye7sHvbiLxlmkfyz73Y3boDrf9QmjYve1ub/412tG3bhtYd+jzb0LavHyPfprRdeGnZm+MYsO/GNm35O361Ddte78KpT0JIyNcbk0fq6iAOv669v7SctnYcfmcc0Yo/HC1X8s8HcfD3eunFwMejGP1sEKeK40YZ/7ucMRvBYHHaKQyWzTP3uqaQzyJ1M4jRdw6jva137tZ3/n4Qg/o28A+t89vAZ1FklwyUV/C7a9t45M8DOLyvDdvkNlt837Y2tL8xgNGb3GBeqD29cJ2Qw0VZjHwdqUppUeZuGOPahZ5yzgl7rTuJd6MI3JXDRW6ov1tuVG6CZV2D+STiV8q2s7YdNYMRa98wnMVnS0nEijkxDIVUAuHikA2Db6qy+ftq7toYTL8uewyX9SNe8XekdSPr3ryEEsK3X1Yek2VVldgWyIjJtxY2hTRaLuiV02SrhLlp0pOY8PXVet98sZ2v31RR5KZFoNSSodSEvUaTzGJZpEJx7l6p9ZhZ2E8PCd8Fn/CUfy5FXUargZdMAxWKc6kpMRUcqqzIrhXr6QkxGZ0S06UK7bMZkYhOionTRqsLvQKn61Ko+HriYnUzXL3o65oV05ddFc16K0qntr3VaZq6ot89HRDayaL4urnPO996LDwmXLKVzNruW9SI2ehQ1fZRr2LxrJg6p/9OVuFLyklVMkFH2XK0UqsV0oqtUYXiYgX9xpahNyJR9W12v0dMpnMil5kS3mLDELNwfVNW8bdYIdhW9+/SEG3/KD8eqKxY/EK8xMFNVcshrTi/WdsD8GwuI6bjATFU3GlKxS28n1Y2o8TR+WaCxYOI4hRj0WmRyWk7XTomAueMVgbzy6h9kJpN+oRdnmRs5yZF5hf5gkY/aRmtuspKneCmvPVYZcuGUi4I+X5lSEwxvlmo0dZSmmLejLJ5a+Xk0JVOLraLZaeUpzkxk5zUguHydVmF/YitmE/EowUmvgtDwt2zMLhVFrQGWenvrl0klLarQ2NiWk6dZ6QfYHCzERZewNVqPSRmp4RHf32R5t+x85XLWXaroUWtTXAz+4NnecvIaMfWC+5ifiD9ItB1ZkJMVbzFCPpq7SvLUn08WNAUntbDyxvcLDgBrV8ytQUHhk63CCQzIqcFPwn9CrvsQKGfxBbmeNCbaJadXLRSkZ9BV3b1XDfoyE0Kd9kyagY3T7R55vJkuBck/ao+gCz4HLSs4Ka8qWix1DzBzIiJHv11lwjVSFFQfVWt39mr/FVmxfSXMu/HXKkKkFf6u5dfLdc8aBsnLs8PPJxvhAXNwmvcccmFXcXXXOH6+/JWCG6q94NVB9RrcddGV308WNO/HdXz8ta5Kdbcr5R69CLqBTgQuDoG534FJpMC64kJTH8137xQORZC4kx1zZ8WqG+45LAh86y85k0B0S89862+jqvoqvVc27RjvpJcHdkbfoyXlnPICkvVQ+uW1zu0bzAvdJ/9U61KWVPRoutBRB/I4ZIHUYQXqw9RwYGBE7aqugYtsLw1itHy9SCC6FyT2FX87j9nETWGgM+0ddzMVtUJU+C8nIO3ZtNiWm/mIy5UHDkWVCyWFYmVITh76tcbXOB5WcuipqQdUy8PI3G6qq6NXrfsSln9uJMjCN5+EecNWq6XuEJxC1pLEYWUf7b07lpsBbJoWaq6bwfMVU3NWxo47hdM5oqTS9v28jclEf3j/A6mvm5evC+YuvJI3DKq1BVtr/F9q/5G2eeNVG+m+hTYj5WffqIY/7ayuWjy23EtFLGWVXJcCQXq0fItSPtt5yoWr+J336kt1xjSRDCstqFVr0Ss6hWg/Qjf1YKd7Sa0vCJnoRdrpx3Oc+UHuqqKxbIisfWDXqiLbFzFdBflvp2ubIzRbB4EMPpHMwbfts9fKDwM49TvduPwyRiU9ycQ+mYCHmUS/bY2dH3CTjI3m5c4uDHDelQOStFkdQeZ1eIY1VuB1C0WdH0YXZOdPv9jHOErwxh48yDa27ahdd8pWZPf0PbrsqussiaMOtP2Vjm0XAXky5tGfnsK7dXfsfxzKBaoppWfbslg6nFWJF2raC5aiCJ0UQt2etxwrLC1RolinktRVmUVv/uew3Afr7pK0PaA1E0/Rs+eQq9VC3Ze78Xo33jo3xgtUN8crGwFesWPUPHuYKnptB3uo4u3E23bVf16HNPVdxibhn4nvPquTRKjx3rh17/zGS+Ge8wwaRec9g8Gi3fG4h874L7GOzibyUsc3JjQcaDyShZfRBBb9Bhsw5BRTwniaQJjPXKypH4Zw/Tn9rlHTCuR/9so+ve1YYelC70ng0jvdMB7fQa5nwIVd24WE3u4VE+7jXJhLDqFqfKS1D5LLofcU+1vkJnG1LsraTxPFVpUHH6nbKvJ+hGNG9FN/noQI1kFnvf7l5lTZCkKOl6vt6Uu53dX4Pg8isD7av3t/kEYgwfsGGUT2I2xXzvxVhyrIsbdQT0jsX7H94QLvXvkS3W0/E6FWw4bFt5h3Gj1g/dl0u/afNYB7wfzd20KNwIYvG0MO6zt83dQd5phKQZAembkgBYCNWAfs9m/CC91nhvlkAvuiq1sXNthG6xDst2KDpscllZ+x8RQuD0C+4FBBO9rBxy9G4j0DCY/16/YtauEZdwgyT7KrNEt0lm0WVSo3WVF/ywmk/Zd5Sy0JmxHyq+uS48OZH2I/YPoP7TGd8gUF2x75fACy/zdt1vgvDiFzNMMpuMhTJxxQd1bffiOI3Bzc50MXx5m9L5VWWdPvzsY/F7PSKxg6HjZo5d6TCocFY+39GX4ESm/47fRdmnboRzUNVLNYCF51+aMB/1lAV/m4fy2a3mt/O/QhrbX5eDdGFK1ulf4OaPtyfOcv1nbyxSq7aUObvTn0Z7zFXk8EflwdIN22DyiX5X1cXXcBUej3UDsVCqv6r8II7qi76AtpyJ7aBCxe6xT80JUX11fCSByI4rg9frp8Jcr/dN8Ckfb2fI6Fqv43fNx+D+LGAfv7QosnQ64Lkxg6l4GuXsTMlmaIVlRCZ5epOpHn8iOoP8PRuBcux+pai1Q3xmv+D2RHYfrg2DFiXtDKVaoZftQNLuCx0SluzZvV9Zvyz/LyaEVKBTK7ujYoVorg0RaHy93cKMxvzWGQHmdgeIO60fqhaeOTyH5hRysRdu5yu/GJNJlO25LB7oqspH64b0U165BKuVvTiIkh+ux/K7y5vPoX8Ob5+DV1KqvrsPw/GEY0brp8JcrjehVWWum04vx05Xh0op/92dpRM76Ea9xxWra68LAW3JE49jFg/qGqdHflG5ZgfOrDngvVnYcmf26H/0fRtfoTvFqmaH2lV2s3ksvs/5j7bs2uvIuFCpb1c5itvTllQ6Ya/RllU2X5YXvcUBd4hEgrY2XPrjRdwjn5xH4jszvstmvT0F9YxiRFQT+K9eClvIKo1/7EbibR+FZFqnro+g9MjDf5HYBE+zHPBUHneQnXWh/cxThm1FEb0YQ/KQXHW8HKw9CT2cXBECmQy54O+WI7lo/Dp7Ugr0FR68CsjeDiLKX2zVjOmCvaLab1a48G2v+XS6Msc8XnmzSXw8X+1FTur2IXR+CtepifXW/ex65mhcDeWTm5rPDbuPt+I20oL+pFQTO5mMBxIKuijvF8T8ehEUdRFA7Xi1KC4Sjfx5A72cLL7zWirnPA0/pQHgtgenlrOjuOIZr3LXRmQ71wyuPz7EHZWF/XrsovWEM2s7aUVVTQVPAdLJUFX896s5RXUa6GxIiJ2IXHFUp6/UU9D4RiCTEjJ4tuFhmRCIaEhNnnMIyl/TMmNf1XY3UTLMzYqIiQZUqfLdyNZKdlVKgl88rS6dbhL4bq0gtbv1oSuTKMhDribBC5Vlkq4pyxCcSuaqszId8IlYjKdxcavLyebXvZzvhMdLxn1CN784MxZWeGttH7FJ1wjxFOC/FjG2oTrcHhuptoH46/HLVycuKZY9DZigudZ+gCPWjyqzV1Vb0u2vvKc6n2IU3Mj2/Tf6SE4m5xIHa9w/Wy31LL05OhE6UfletrCJTbu7OmHDsKVuWLMpeu3B/5BMT4VIXHBPatuMRrrkuYMqS4v0yWzymzsTH5hOQlpZzfEzE9G4Rcsv/hDNflhIXKmIo2uj7c2LytCKUc/X/JnP7h+IUgXRxikhcNNalHJ/PMl9hdkoMlb7bIhmgae0xuKmmp+Q+Xx241Cl6303agX9CP6gv2CMWZt1cUKozBKdDwlV1wDCfmBDT+gmxdBKpKFWZb3+ZEaFSv1JlxfbupNypFnY5USw1MmbOpgLCUyNtv1EUYTnmFaF7zE5cbkEW11plkf68ipI+YS3N22Am05rBTXnRAh1ftLGcqMv+3bWTnK3GSa5U9G4g9GCeNof5/qaWc+KvRwtMImM1u/dYULRjpfN8QMTKN8Oax7TqUtbvXsNKGb21ciJUla27Dm2/s9Xte6tMZkqMnbbPnR/MnS7hDWvHf/lytVL25+KxuhgQ0YuyTf9H++NTDYVsCslUBrOFLFL3M8CrFlj0Lm1b29D+GzOU9cjx8jyP1N9jSN7JFJtfqp3mpVsyVPs5hfi9FKbvF9B2QIV9f+k+7fIVtGUl72WKrQUyLdr3N5vRvtcCha2l1kkSo9YODN5VMBRNN5TZN3u1F23O0q1vBwI/+WF7lkAqlUZe6YL9gAWmZSbRW+7vrvd8n0qmkdN7wH/cCsteBW3mLlh3rcM+QptPIY/sQ+2Yk5ktHjf1VkNt+vbS0oo2iwVmxbTgUc+6+jmK4SMHMXLbCm88hqHOJdb+LI/88xaY1vKY/iyO4QNdGLlrg/dWBEP/uuwjOa0CgxuiTSWP8Mkd6E36kEh6GqrsuSC4eRSCs0bFRqKXSinAgRexG0OwvdALsgLin6jo+rwVQ1+F4D3EwOZFY4Vios3kQQj+K2vX/JvopbVThffGNCYsAfivv9DWIdoVRwT+b9oxcXOSgc0G4Z0bok0k+VkHOi52YfLeWMOtpHjnhoioEu/cEG2QwoNUZS6Zh354zyZhO+taZvNvIiIqx+CGaIMkv2pHlzoI//UoIlcGcdB2CuEaCfaWMjvXuzcREen4WIpoQ6Th//1unJIJwIr0/sTiATgb7HYj/2MUiXQa4fOnMC479dPZ3p3A8FEzWto6oP6Gz/uJ6OXD4IZoQ2hByTuD8N0II/7MAvtRD7yfumBdRiyS/PNBDH5vhvWQBW1yWtHTNGK3U8gf8bHXdiJ6KTG4ISIioqbCOjdERETUVBjcEBERUVNhcENERERNhcENERERNRUGN0RERNRUGNwQERFRU1n3puC5XA6PHj2SYy8ntrYnIqKXiaIo+Jd/+Rc59uKte3ATiURw/vx5OUZERETN7tSpU8WyUdY9uInFYrhy5YocIyIiomZ35MgRvPHGG3LsxXshGYpfwCqIiIhok9i2bZsc2hjsfoGIiIiaCltLERERUVNhcENERERNhcENERERNRUGN0RERNRUGNwQERFRU2FwQ0RERE2FwQ0RERE1FQY3RERE1FQY3BAREVFTYXBDRERETYXBDRERETUVBjdEtMkVkH+YRPjKIEav5+U0IqL62HHmZlbIIvrnYaS6J+D+rZxGGyMbweAfRpE0WWD7nRk7nqYR+9GMwcse2LbLeRaVxPjvBxGGCZbOLph/pU16kkLkbhrY60Hgoh2KMSOV3Pej9wM/kjfi0P5KRY5gBqFj/Eu9FPJx+P8aQ0Ph7F4HPD1mOVLl5xTit2KIPZBLammDtccOdY/JGKfmpAc3tLnMPoqJwHmH0HZVPfAU3rh8gTbQrMjlciKXjglfn/G76EU5HhAzco6lzOrvz0yL0Dlb8b22cyExndGm5WblHFThF/k3T00Ih/x7a8GNfJGa3qPA3O++VKm9XcyKxCXtOKqownUpJKaiU2IqPCZc3Yr2HkWo7082vO/S1sPHUpvF8zxS18cxYN+N1t+5MJkqYFa+RJtBC0wmE0y7bOg/qh1ypezX/ej/JI6CHF9Mi/5+xQLHSRdUODDwtgMWRZtmapFzUIVX5N/8N+3okJPo5WPuVKEeUrV9RU6ooVCoPloWEP9ERcenLfDGpzDxrgNqt7aco25M/BDDRF8W0T8dRv8fkw3tu7T1MLjZFLIIn9yB9vciaOkLIJOeRuBCP7rkq7QJ7bfBJg+28Y8dGPx+GXVBWlpRvCH+SnGMiGp51VwMal1/msLUD1OYzgj9SUPNMnmi6pHU3XF4Po7D+sEgnLvktJJXzHB96INVG4x/OIzAA2MyNRcGN5uCAvVCBrM/TcJ3wgaFF/Kb3+sD8H/rha04ksX4224EHxZHiGjNqFB2ysGG5RH5chBaaIP+Q3oIU8N+Ff379YEIxr9PFidRc2Fws0mYFAWMabaWHZ1DCASdRkXgbBD9x0YQf1Z8iegllUL02xTyz+XoqpnQutwDYz6GyBf6gB1dxQCmFiu6jhhDyVupuQrr1DwY3BCtgvnYKMaPy+dTt4fheCeMrDFG9BLKI+5oxw7LYQxeiSO7ZkHOMtyPYVz//5B50RaIZousO/dtAmlWvGk6DG6IVkWB4/MwvJ3GWPZrNzxX1+A68HkW8asjOGXbjW3btsmyG11vDmD8ZrZ+JchCHum7YfjP9qP9D0Ej0HqWgv/tdrTpy2g7hfBc9DWfP6Z/Xz+Cj41pyT9r722bX+fhs0GkyqsUacsLnj2M3fJzte3TTmRXU0tWzMzeDmLkZFexlN5b/E4nRxD+kWeX5mCFI+yFsyWBUe13bnutHf2fRZBeSXqibAbTcrDoeQHZH+OI3owi/qO2D9QJnNLplDFgakWrMVSTaWepnk4a2Z/lIDUNBjdEq7XdhqGrATiN51MIOvsxcnvlJ+vCj36csrTBEzfD/f3MXKXJXNKLjkfjGFDb0H7Sj1TFI7AkRve1YVvrDuy29uLUZ1pAUvwIaQTfUXHqrykj0Mm2oKVFm/tPerDTih3mDvSeHEXwvjbzc+2qu9jCJKkFLCrUTv3gn0bkMy3Y6RlBUl9ePo6RN1R4rhdgPqTCtkdb5P0IRp3tUOu1PMknMf7mbji+nEaLpbdY3Bd88JywwawtP35lGL3dLoSLwRVtbS2wHB1C4F4amegYXJacEQjv2I3eT4KIL+e2phbMFGCB0pJC+BNtG3ytFW2WLhxUD6JL2z9a/0Fb5mfRBXeHsg/DxsBO06KP+lt+tUMOaZ+VwU3z0Q6atBmV5XhgnpvNJRN0CPQFRHVmjdx3bqHFN0buDcUpAmn5QrXib+sQgUdyvFw6ILQgqbj8mjk4Sq9r66idYycjQifkZ+gbE4HzqlDf8olAcEJ4+szadK+IyTl1mW9cxrywCvsRu3BdTojcL/JFzcxXzrnvZD/nE+5uu/BGy775LzMicFzPG6LPYxcTP8npZYp/r56JJb+P9WJCTqwWE175GZjnZuvJxAPCW9z29N9Q0bbHMTH1qIHcTsX9RBGKtn0oe7Xt+IRHDJ12CHVvaXszSuV+MCumzsjXzpdv6TXEvXPL4DG2+TC42awY3Gxa9YIb/cAaO28k6CuWegFK3eAmIwIyQaDnh/oH/1y4FJBAuMI5OXVe7Lxcv3ZicAbLP4EWiJwYExUhRNkB3vFVrU+bEL79peXVDl7ETxPCLpfh/GbhXyUTdIqhaP3vM/25aiz/aL2EiAxumkHup0nhO2aZC5bNfV4RSi7cfucU9xObcAenRfVcuXsTwrXHWI5e7F+Wtpz5fWjJ4IbH2KbGx1JEa6YFto8CCJQqGF9rPMFf0d0AfNf0ASe69tW/oW466oJXrsL/RQh1a/goAxg4Vp7/wwznZXcxv8dCDvR2V+UKKbKiq5Sz8Gg/1D1yuNyeDqiyVUr0QcYYKKMcC8DbXf/7WKyqMfBtmpWxm5hpjx2e4DQymRgmzmjh8LVh9Fp3YLe9Tj2y7RYM3otg7JjFyAtVxrTXhYnrE9CWUhT5OIi4HF4+bWdizqmmw+CGaE1pAcRFP9wy+NAT/LkarGCcvhOBkXGjHW3y/bVZYTsuB29EkawXERwwa59m9ZRdMrrRTgC1K2i2oe11Yyj7vNFQroB8Po3kzTD8X8XkNHopKDa4Lkxi5kkCofNOtCSNemTmfaNy+5dM2na+d5H+n/b0w3NO7ijZBNLVdbbuLREsz9XV6YJ50f2NtiIGN0RrbacdvrIEf0HnAPwNZEHNZqNyaCkt2omgdDulULfVyFppXYur2lL3Im92Fcvutg70fjCO6EMtHNzXJmeil8pOKxynB+E9Zty5y97XKxAvRwu6uvvlcFgGN4q2PRUnLC2b1t6la0EL79w0HQY3ROugpTzBHyI4dXxtE/ytScDxghTujqPXoncvkoTlg3CxTGemMXVZbzHlgGpdi/tLtJUUHkcx/s5h7P61FuT+aRrqW2OY/MktLwgaN9/iyQHzq8ZQyyvyNkx+trH++ZR23rlpQgxuiNaJ+dg4/KflUbOY4C9Yv35MhQxyDecFaYNpuxzcjB4H0W8dQPiBvVhHwt2pFMtiTXSpeeXvhjHy5m60vnYQA1+kYDkzgVgmg6kv3bDvWeQR1FKUjrngpm2X7JXvxuKPpbKPZBado9Y1eXxLmwuDG6J1Y4L9QnmCv/5FeyGeq9uCEJI/ysE6Zp/J6KfHCssqzgnrLXnNZ9z6r1cZmV4CBWRvjuOU2oYd1l4M32qBU69zk5vB5AXXXAe08/JI3W/sMiD/xKjAbv3APnfXR7F0yUrzcUwv8jg489hI9uc40LGgwjJtfQxuiFbiuRZgyMFFVST403shdmP0Vs4YqWLeVzooZxG6sVgrK+3gn9Tr5yhwn+7d1FedhWeyimjdyshA+qeEHKKmIrNs977eijZ1AP5HHfBcjiGTnkbgjB3muhFFDrGzgzJj9mIKiN0Y13YDD3yny9oA7u9Cb3F/iyJyp969myxScX3bdKD3AJ9JNSMGN0QrsZxmy7uc8H/lkfVv4hh2DsiKjFX2u+CVj7GSH/sQrtfL+OMI/F9ox/Tjo/AcWf9rzrm7RKtxK4ZEjSyw6W8HMXxlYfNx2sLyaUT+fAoHX2tDl3MYsRYnvOEEcqlJ+E7YoCz5TNKM3rd2wHMpumgF48LtUXivqfB+64Va8WjWBvtZ4z5O+OvJ2o+CH0wicE3bh8654ZCPs6i5MLjZpPL3Eig1kA1Fotq1Om0WmYf64TKNTMPRDdByyIvw+aWqS5pg/3AcruLjm7DRy/iCHz6N4IcexLq15X3urHHXJq9dGZcGc9o18OIKuVJgkUeuToXn7CPZiutWeuk6Q09yFSck678OGUFdVvtebw7Cfz1aLNHrQYw4D8PfMojAp4eL8wI1mvPqfs7MrTevfSfazOIY2bEbh9/zY9riwlhU++3uBTB01ArTMirBmw7Ycfizfrj+mlzYw/jzPJJXB3H4rQT6r09iqHNhtGR9Swt89MfB14cx/HX1VqvtQ+eHEenU9qEPVdb/alYymR9tuIyYfF8Vavd8Bs/qUkxB/v5kjcy4tO4yk8JzyCa0YGLhb3KpXrcB1cq7KqjT/YLuSUz4SunqFVW4LgTEZHRKTIXHhEvbPhznJ8VMWRcJhoQYO6QKi+zKYL4owtK98DMmLmnTqtLY68XcqQrPd/oWZizPVpYF1ihmYdOmj90xlmMoywqrF8Ui1EOlTMizInHRvmCbVo54RexJcYaKLMnYY5t/752x2vuDPg/3g01K33aXyDzckJyYPF36zY1tTi0WfR80C/uZgJheahU5bRsu7kfaPnDMIybCU2Iy6NP2IUWY+3zz2x81pW36P9oGREQvwvMC8o9TSCQzaO22w7bIU6W83gNyPIZ08XGOCeYDXVD3W2Daipeaz7JI3Z9GpqW9ONq1ny2maCkFZO/GEEsmV7UPVOxHO83osqmw/YZViJsdgxsiIiJqKqxzQ0RERE2FwQ0RERE1FQY3RERE1FQY3BAREVFTYXBDRERETYXBDRERETUVBjdERETUVBjcEBERUVNhcENERERNhcENERERNRUGN0RERNRUGNwQERFRU2FwQ0QbS+8x/GYQI3+KIC8nERGtBnsF32QKj5OI3IjKLv51JpgPHYZ9v4Jl9PJPay6J8d8PIizH8KoVdkubHNHlkL4TR6rs7Gzeb4fltTZYD9nR9RsTf78KKfjfHIA/GUX8gZzUF0DmGycUOUobZ8FxaKcZXTYVNm07rif9/SjCP8qRRbVBPe2EdbscXY18HP6/xhoLivc64Okxy5EqhSySN6OI3c+gUJxgHHcd2nGXtig9uKFN4Om0mDhh1gPN2mWPS0ykZuXMtBFmczmRS8fE2HFl7nexn/YJ36WAmIxOiSlZQpe1aR+5haPbIrRDo/z97MITnBY5uSySf8+ctt33yb9RX0Bk5Gu0QZY4Dpn7xkTiqZy3Sux87fcsLF4Rk+9ZtUcB4ai5joXFEay9dc185xGqYhaO83I/jgSEt8/4Gyz2fWlzY3CzKUyLsU65E+6xCcfpIeG7MCTcfTahXWfM76CKUwTS8i20YTLfOOd+E29cTqxnNiNil/WDpzG/0j0kpp7I16ho7qS4BsFN7t6kmPpJjtAyzYiJHuO3MPe4xVhYBuzhMeHqng/o0TOhzblQ8XdULEI9pAq1c5ELNWj7wC/yTaslgxtzp7ZObb0WuZ/VKvbLCz/1TNCpXYDYtP24+sIxJ6Y+shbfpxwPMejeghjcbAaZkHBqO5HtfGzhlb1+ZVt2pwAnQrz632hx79zvsWRwU5IOCGfpwNupXbnyanDOmgQ3P4WEu8c4oTb8m1Cl3KRwayf6oUiNX+GXGRGYOw4pYii68C5yJujYgLtvMeFd6W/+RPu++j5Z75j6JCRci3xf2txYoXgzeJZDdr8Xo2dsWPBE22SB6/MwhvbL8SshRLNymLaOXU74rw3Bqg/fHobjnTD4M66hn1MYv57WBlQoO41JtEw/JhE648VwT416Jq+Y4fzYC3txJItUtk4tl13KwmPYulvZb578ehjj2k7oPqrW/sw7VThO6wNZjFyLsrL7FsPgZhPI/wy4rnhgq1fjdLsNA2cdciSL3DM5SFtKy796EQ1r14Ka7Ne9sH4YlZUXae2Y0Mqa2yvTOYTMBbV+xfc97egqDthh/12dira/at2AivMr+M0fBuH9IKkNeOA4VC8cM8H+B68x+MUoQqWK77QlMLjZBEydLjj3N7p3WmF+VQ7SlmPaZ9OuMw3Z60mk5LCukE0henUUA/bdGLltTMvfGMHh17dh27Zt2P2ZfjCu8jyL+NURnLLtLs5jlN3oenMA4zezdYOnWutaoNRE+2QXdn8SlxNrySN9fRwDb2rzlT5DWzsOvzO+wruMecQ/60d7W+Wy4nMtCGlD3I0hpP1n+9SL/j3GpHKZhzE5ZCjk08UWSNGbSaTz6xTGZzOYloNFzwvI/hjX1hlF/Edt+38up1d7PC1bPu5YPDDaqcj9NYost7+tRT6eok1tVkydMeolKOemtDHaUCupc1MyOyU88r3aVaOY0n7MzHeuyorjWtGXO6utx1Y2zXoxIRdimE1NCNceCNu7AZEoq6ScSwaEW1ZQN5+YENNl9XvqratCJiRc1RUzz9dp3/JkSgx1K8J2blJkShtmbloE3rXJ95qFTa9gWix6BXl7sUyUVYyvqHOjt9Y5rghlr3zP3rL6ZjUr1OfETMQnnOXzVZShtWuZ8zKTdW5s2vGnXp0//Xd0fDUjMtGxufpP5UXp9ojAvTWuMVisUKxX0p8WofPOGhWKzcJxYUpkqiowF+sH6a/v94nKvapKWWss93es7biV8M7NVvAshsjX+oAd3pOL3Damza+lVbtWLEkXrwaVIxOY0Sv33/EadXI0sz8F4Xorga7LAQQuuKAq2mFaKcur81B7vfsU/Fbt9UtOWMvqHJj2OzF2NQCn9p70lVNQ3wlqazLUWtcCigMTGW2eX2IYqvP0oehZHCNHDmIkN4jxT+1QShumyQLnpQAmevSRNOKPTMVcP/ZDvXAHvcXSv6s4Z6X8JIado8AHKWTuTWHqB63cy2Am7IYW5GlX6UF4vqx8lJe/HUL4PtBxoENOscJ5zgffBVkuH4ZFvkIrUUD+xwhGHF0Y2xVG5NM69VOk2NkutKmjSDy3wv2RB65DNpQyy2RvjqJ/nx0jt9f6Lo5fW2473N9k0dbjwdBpB7SgWL6WRvjsQVhPzu8DutlnsgbN620oz1a1mMwzPkTeUmSQQ5tY4qJ+FawIZ7BWA0x64VZz50a27jDeb9y5mVOes6OqRdVs1Cs835XaoWREQOaG8fxQ/z5eLuyS64FwhauuOsvWVf87zK+n1p2b6UvG3Rnl09rXvvPrd4vJRS565+7cwFEn1YF+57J0Z6ZOjpS530RbxiM5jVZs9geP/HtWlT124b608E6ITv8dlSNeMVXdXOqXjJj8SJ3P+aRo2/1atRYsbsc24a6RQyp3z7izWfrs9i/nj5+Nt9Cb31/r5cmhzYl3bja7u6NwfxCH7XwY/mN1smvS1rS/DTvq3IZzvOeCrSyDa0v3EHxH5NXo3QB81/QBJ7r21b+PZzrqgle+xf9FqOLKdfWySN4y6uF07ap97WsyW+TdoQzyjVSC7+uFWuuODlqgvmFUxAYSSD+Wg7Ru8q90YDKVgxY66xfAEE+mEfpIhfIggvH3DsL6+1EkK25kFGA+FEMsPFS8y1jhFQX285MInJEvZEcx/v0atRXcbsHgvQjGjlkW3FEy7XVh4vqEbOEFRD4OYrGaY0tpeUUO0JbA4GYz02/7nx5E+ngAgY9sfBzVDB6ntdOzdMha95FJh7n6DDEvfScCo2pxO9rqz6axwnZcDt6IIrlObc/T2YwcqrLdJB/BmVffPPuVVjlAL4LS7YS9vMuQnRY4tAAl/JERrmZvDmL4anm43AKl0wZz3QBAC1DfLjUlB8IPtPfKLg/0yr/LLn+XleVN2ja+d5EHZXv64TlXCqpqBMbX0ksG/bPy//bXFt3ZaJNhcLNppRF8x4Gx1ycQ/dw599yatriHpVYagPtAx4oC1mw2KoeW0gLzvlIKgUL9liMr0gKTDFiSN7SThjFYIX8vDv2T2j7tRxcj8ybQAtsf3PN3Qm4ml5eraY8KxyE5fE+vbxaFVz2IgyspFxrNO9OCru5+ORyeC27Me0r7xRK0ixGjRaMW2PDOzZbC4GZTyiP6oXbF8XgA4c9dsKxFB3O0KSTjEWNAGVokv8baaV23A7IJat+Q0cnl9WEMf10V3jyLY/yCH8rxCfjfta4oiKNNaC7Xjaaw3Aq2rWgtbfL79Lt5KoajU5haSTm7eMXmci2/KlXhd8yn0XiltEXmMNvQ11BhYQqOLYXBzaZTQPwTOw7ePIzwd0MV9S5oi3sWReCi8UDJdrYX6qrP+BnkGk6b2gbTGm9LLd1exMIeqEoWwT90od05iNHPRjH68QAO21yYdsSQ+orBedNaRTZixx6ztgEpsHarUFdSfqssP2BWOuaCG+W1dmMAaWQXu/2UTaOYvWe/9l4+ldpSGNxsKlpg88fDcEQOI3adgU1z0X7bC9rJXz+Qdnox+lbdhthLUnaVbqmHkPxRDtYx1+S1xwrLss9Es5hdIngyH/VhKp1B6P02WG1dsOy1wmp3I5CcRqBWdyK0KRVuDGPkb43diTHqoCgY6ukyAoxsEqmGEtxlkPlJ/9+F/kOrjRTySN1fqraMIf/EqBNm/cBupBTQ/cYKd3EgjMRP9b93PpMpPnpTHF31UyfQpsTgZhNJf+2C47IZ41e1wIZnhSZSQOqv2m/7SVI7SjoR0H/fVQSu5n2lA20WoRtxben1aCeApF7rRYH7dG9lva2dytx45kntCKZwOwj/DTlSz7MU/O90IWANIfCuA/Ye7aq60wIT6ydsKS1mM2KfBpZuUXc3huKD1b5RuLqNeyeFVAiuLxbbDqUftffeBexfDsGx2grmyCF2dhDBJVvOFRC7Ma7tAh74TpeFJ6YOdPUZg4HbNTJ/S6m7ek5mBQOH5sIi2iIY3GwS6av96PpDEr0fOmFK12gdoJXI1XHjtr9WIuznZMMUnubkUAPyKQTfU9H+dhBZxQ7fdT+cNZs7L8N+F7ynjSvf5Mc+hB8WBxd6HIH/C+3QfHwUniNV0XLLDuyQF8/jX0cWntQeBjFwPg1zqQJoTWn432zHqStpWF7dLFXe8+x7bSX2dEDNDFe1gKqWRvCiD8lONyYvzDdyaOnuRW/Yi8CixyTtvZ96kT4Rwthba7GtmNH71g54Li3eP1vh9ii811R4v/VCrbigUGA/7i7WGct+HkK01jZTiGLy8yzQ44WzU06jrcNId0MbaeY7d0Wa/UbK8pPH0dqYFbGPrHO/g+vzKRFLZUQul5MlI6bjU2IqPCGGTutdDejzKUJ9PyCml8renvQJ7dqyuFx3pH5yvqJHofkEZXrCvwXLNtLlK92VyQDnzYqpc/NdFijdLjF0OTT3uS1HfCLxtCyJ37uTC7v9KE86WOxmwSU8F3zCV6NMROslQJsVk+/KZSyWCn8uSZ9V+JJyWrmKBIjad9HWqX8P+8XEws9NNcU+0v9+tRPiFRPxnbNp24ne1YGcViZxQdsnOofE5KOFf+3ZR1PC12fR9oFJbatcQ0+0fUBPbvplQuSqkwr+khOJoEeoex1iLFlvCzD2EX2bsZ2PVX1nbT8/b6vT5QdtBQxuNlomJJylg/IyCoObjTAtxrrr9WFUVRTtYN7n1gKGSTFd42RQLvOddhDuXNgXD/bYhHrIIybrxQVPYtpJQ75PUYXrQkBMRvXAaky4ui3CcV47mdTIJDvnaUKMld4/V/RALCTflxDe8r569M+jnaDmP86sSFxyLOirql7RTyD6acY41WjrPqTW7AvIVus7l2WF1ou5Uy3L2KzTTkaflmXB1YpyPLC2J9Mml/nGOf/3035rx+mhYmDqOaH9TnvtwlMr6CkpC8zn+gWTfYMpe53CVze4XY2cmDxd+r2r+zAzC/uZBi4ofpkRofeN7cbc4xa+4GQxwPdo+0WxL6wUQ+Otapv+j7ZxENEWldd7QY7HkC5W6jTBfKAL6n4LTA02J9F7b04l08i1tqF9rwXKcuoD6Y/dPnTBm++C61AXrLvmH3/pPY/HbkUQ/KL02MsKX9JIYejZX/xvzc19lx1mdO01M6vschXySN1NIJVKGtvTTu3vaFNhLU/oV8/PKcTvJJG4nyk+Kmp51WK8d9ey2zUtQwHZuzHEkvLzrmD71xUeJxG9FUPqsfbJW9pgOaBqy1hBiyzaNBjcENEKFJC8MoDBDyehXIwhsFjXIA/96DWfKiYvdASNliuhY7LCDxHROmCFYiJavtujsJ/0I2r1wrtUn2e7DqNXtkwhInoRGNwQ0bJl0wkj9b6tfRldgyjoeN0oRETricENEa3cL/L/xRSmMX1L+79zEPbfoliIiNYTgxsiWjbld3ajA8W/1skRUiZ9bRwjWRu8f3IXkw+WpVIjIloXDG6IaPn29MN7yQ4lO4r+d4JI1+xxPI/UX/vRdTYP760IhjrZ9oSIXgy2liKiFcveHIXnPR+ChQ64jzmgdpthKmSRiscQuRoBjnrh+9C5gn6tiIhWjsENEa1a4WESsXRZtxTMM0NEG4jBDRERETUV1rkhIiKipsLghoiIiJoKgxsiIiJqKgxuiIiIqKmwQjHRahXyyOvdIEstpgZ6UCYionXDOzdEK5JH+vo4Buy7sa11B3bsmC+t29rQ/sYARq9GkfpZzk5ERC8Mgxui5XqWxOgbFuw+GUZLXwAzTwX0G6B6mc3NIHa5H213xjHoPIj2C1GU3dQhIqIXgI+liJajkMRIdweGbzswkQ7BtUtOr5aPYNByGNGzCSTeZ29KREQvEu/cEC1D8gu3FtgA1gvD9QMbnUmF/ThgVtrkBCIielEY3BA1LI7IhXhxqPfAUndjWrBjh4K27axaTET0ojG4IWrUwzQSWX3AAbNSnLIo6/sp+HrYYyQR0YvG4IaoUdk0wsWBMBI/NVBNuMXEjiOJiDYAgxuiRilmOOTg6MejiD+TIyuUvR3EyMku7N62Ddv00taOw++MI1q8O1SSxPjvD+KgVk6dHcXoZ6MYPGmMH/z9uPbqvOSf5fSTgxg9e6rmPCWNrXteIZ9G8ttxDLzRjv6rcqaHYQzqTeH1978+2vB6lloXEdGq6a2liKgR02KsE3rrQqPscQhvcEpMZ2bl6w16EhO+PrNQjnjF1KP59+biPmFX9GXbhPu7GTlV88usyPwwJKyl9R4aE4lcTtRa62xuRgROGMsYisyIXPVMy1x35juXMJfWK4sjmBEiHRDO4vyl4hUx+Z6iRdZT93sSEa0RBjdEyzAb9wpb2Yl+vijC0q0K9YRH+C6HRCxVO/gQ6ZBw7dHm79SCgadyWplc2CWX5xCBR3Ji0YyY6JHr2u8TCTl1oVkxdQZCOTO1cP0rXrfm3phQ5Xd1fB4QXu27ui4ERCjoEy494OsLCC3kMSyxHt2i6yIiWiUGN0TLlZkSvmMWociTfd2yxyF8t3LyTbqMCPQZr7nC5dPLzE4Jj3y/9WJlCDMfEFiFLyknVtPeP6TYxcRPcnzO6tYtREx45Wt6QDKRlpN1T6eE9/1JGdw0sB7dousiIlod1rkhWi5FhSc4jcxsDjPJKYQu++A744LaaZYzSA/CGDxgx8htWfn4bgC+a/qAA6q1TiuqlnZ09BmDyXgK5dVSTD1ODBVbaSXh+6Z25uP8jTBGelzo3SMnlKxy3eWUc270l+f42a5i6KIdxY/WyHp0Da6LiGglGNwQrVSLCeb9KhwnPPBcmMBUfAZCC3imIz445oKLOIbPB5DWhtJ3IvOVbuu2omqBaaccvJYuvm9Oi4reD4z8Otk/BhFZ0G9VHtFvY/C950B1WLHqdZfp2tdet2PQxtaja2xdREQrweCGaC1pAY+lx4PQnUl4SrlwrqeLdyay2agxjjD6X5OthxaUHTj8hQJLtwr1fSuq7gXBemxQttjyw/99VUjwIAR/1g3HfjleZi3W3YjG1qOX1a+LiKgeBjdE68Fkh/u8KkfyyOXlYJEbk7n5zjYXlgymo1OYKj3qKfeqA+5zxtTIpXBF8+vkt35Y3u9fIlBYxbqXZan16GWt1kVEVInBDdEyFPKN9/Hdul0+HFLa0FbxnCiDfEWwsxwtUN8cRPHh1F0fQjfl5ylEEQr3ov/QUt09rGbdy/Gi1kNEtBCDG6KGZRF+ux/Bx3J0CbPPjf+VEzZY9P93lVIAhhG7t4oz/34H3D36QBYjX0egLyl/PYjUe/1G0FPDmq17CS9qPUREi2Fws9kUssXHC3omWqP4Eb6dQl6eKGmjpZFZUJG3lizi3+udNdjhPakWK+Ca93XNBR/jWlCy8hZCZvS+5TIGr/gRepCE/88K3EfrP9xZu3Uv7kWth4hoMduE/vCbNoECUlcGcPikv3bLkT0OjF0LwL2fvUxvnDT8v9+N8NFpTJ7W78XUl77ajy5nFGowhsCxUi2YLMJ/sKL3a/2Ur8BZ8VqV51mkf1bqd9BZiGLYfBAj2qKsnTa0OMYRO1Pvvo1uteuOY2RbF4a1IUcwg9Cxeh9sGevRLfU9iZ7nkboVQfTvGZn+oAVte7ugdluh1DscPohg9NuUHFlcW7cbzt/yuNp09OCGNt7MZbuRIG2PXbgvhcRUdEorITH2llqWLK5WcjZ6cUqJ7BShvj8hYmXdCsyZzYipT+3ab2YWrsvT9bMEly8nXTbXbE4kwl7h2GOrn6hPSlywyuWojW0Xq1l3OiAcxfdBWC8skXRvqfXolvE96eWV+2FIqBXdfJQVRRWeet13xL2131OjeOPyPdRUGNxsCjkxeRrCdm5SZH6Rk8rMfOWcC3CUczXS6tOLMTstAqftFX0tKXtVoR6SpVvPWqwIyzGfmJrri6CGXEKM1c1wrL/fK0L3FsnuW/LThLDr7+mZEA330LTMdWe+88jvVT2vWdi07+z5rs4XXXQ9elnG96SX0mx0yNh+FItwng+IyeIF36QInHeU7YOK8PxQ44hYDG60bUzvEuWQbUH/aOVlKCrfQ02Fj6U2hThG2kKwpX1Qa94dTcNv341T17XBvgAy3zjZdHYjPS8g+yCJ6dQ0kg9kpdmdZlgtFnTst8DU6B3un1OI30kicV+/3a7fau+A9Xc2WErJ7RpQyOdReMUE03Y5oVFrsO6G1FiPbl3WRU0l+acO2JNuRC+5YKnKSpn/2wjsB4a1I6emZwIzEVdlCoTHQfS+FkLvoxCcr8pp9FJhcLMZPPDj4JdmTF4wKp7WEv9kG7o+1gaOh5D5ysHghoiaWBbBN/sx++kUXNVdiRQVEP3QjIN/1Ot2eRETQ7AZLxiKwU0C7tl6F4zU7NhaajPY48LUIoFNObteiU4OExE1hwKytyNI6rFKcXQaibwTas3ARteC9n1dxqB2QKx97NyBVgY2Ly0GN1tCEjG9VXGnF97FWp4QEW0peaSvj6J/nxltNj9SpZQXLSp8P1Q9aqqjlEeKqByDmy0g/bUPvpYhTH0/BCuvRIhoq3ueRfzKIA6/vgO77YMIFjrgvjQIe8P1Y/KI3dCu+BQnxt+pcdc7m0ZMDhYVtCDqbhTRm1EkH+bReJ5x2qoY3GxG8uql8HMKkT8eRtfnZoSve6GyAiYRbWX5NCJ/PoWDr7Wh6+QoInDAG04gl5rE2Lu2Bb3Z15P/2zh8152YuOmHo1ZApB1Ds0fNULJRLfg5jN2tWhBlPYiD6kF0mHegte0gBq+mitm9qTmxQvFmU4hisPUgRuXoPDPspz3wfOyGyko3RLSV/JxE+AsvBj8OF5OUKt0uDJ8ZQn+PubGARr/ge0Urz9KIXxuF50IOA98H4KxXJ+f2CLYdHdOCmyxaO1VYrTaYC3HEbkURfyDn0djOxxD9yNZQfUfaYvTghjaRzJQIRKZFrpS64ZecmA6XJbJSVOG7w0w3RLT5zT6aqkhEau7xiIn4YkmgapkRE4fm89LMl0VyJel5bhS78EYXrisTKU8MWCdPDm15DG62iNlbQ8Ja2qmXk7SNiOgFyyVDwjuXxFEPQnxiMrXChI2z02IyGBMzT+W4mBWZ+FhZFmyzcIUrg5jZRzERS8uRGmZ/8MwnmOwLiOWGW7T58bHUllGWyA9OhDIBOPh4iog2lSwib3fh8F+LD5+gvjWMoY9dUF9dhwc/D4PoNfdDb0gKZQhTae8yctqUH09r5MmhLY8VircMM9rn9r4CCuwlnIg2oXxeD2wAS9/g+gU2ul0OuM/JK7ysH/H7xmBjzFCPqnI4gfRjOUhNg8HNlmSGwpZTRLTpKHBeiiFw3oHCtUEcfK0Vu98cQfjuerRLKkvkh+xcK9NGtW4vVWXugJldNDQdBjcb7VkUw3+MN5Z34RfjP+WcHV2s3k9Em5Fig/OjEGZyM5i84ETLrWH0WnegTT2F8ZvZdcoxo678gq/P3FCyQNpaGNxstO3ajnXLi0BZ88Takoh9r//vwOjbjXXVQES0YUxm2M8EMP0og9hlDzoe+TGgtqH19V6MXI0jW/dOSwr+D4LFJuNLmS0to0d21fA8i+SPjd0lymSNNbiO2dmlTRNicLPhzOg4lMHw+cV35vRVH3x3bXB/54Nzl5xIRLTZvaLAdsKHyVQO0xHt+NUSw7CzC22vHcSpP0eQXhCLaBd8LcPw31zqHk8eyZt6dWIrhj7sN+6+PJ9G6MQ44kveHkohdj1Z7FF86GijqQNpSym2maKNdWuo2CTR9m5ATNdoLannZbApqhj6YYVNKemFm83NiER4THi+mpZTSPujiJlkqFjGzmjbupxML5tZkYmOCXePea4pt+N8QMTK2mPnvnMLdHpFbK7590Kzca+wQRH2iwltifMSF6zC/uXiyTJmgk6h7HGJ0CLNxWlrY3CzGTwKCedcUimzsPW5xdAFn/CdcQl1r0XY9RMB45otYFpM9KnCMvdbauV8TL72ssqJyXOqth0r83+TYtFOXHIOenkV8+H0lYIchwg8ki+kA8KhTVOOeMVUjSQ0uXsT2jHTLFyXpysCm6KkT1hhE0ORzMLXZjNi6oJDWLo9YpKBTVPbpv+jbVi04QrI/5hEIjWN5AP9Pq0J5gNdUPdbYGIFmy2jkNc75Usj8GYHBm5oE87HID56uTNoGH+TAhIX7Tj4SVJOZW4Rmld4HIX/fARtH/tkX1FJjLzegeFiXUQFlu7DsPdY0PY0jdiNKFJmF8Yueup0RZPEqLUDg3e1QcUCdV+bMfl5BtMpQP1gDKPvq1D07hyoaTG4IVoH8U+2oetjbYDBzZzs1V60OYsp1zRrENxko4g87oD9t6wz0awKD5OIpVJI3c9o4XEL2vZ2wXbACvN2OUM9z/NI/T2G5J0UMnr9m5Y2WA7YoO41o4VBzUuBwQ3ROmBws9CaBTf5OEbfccF3NYWuYAahY2zrQkSV2FqKiLaWZ2nEtMAmqw2ad/KuDREtxOCGiLasHb9ihTQiWojBDRERETUVBjdEK1ZA9uYoTtl2Y9u2bVrZja6Tg/BfTyEzK2epo5BPI/mtH4POdvRf1R+waNN+9OPUvrbistpOhouPXSrlkb4+joE32tFWXJ9e2tCunsLIt0nk62V8fZZF6mYQo+8cxu5P4nJilefad/kxiuAnp9D1+gjqzFWUfxDB+Du92nxln+GNgWJq/eUrIHV1EIfnlqX/DUcQ/nF9kvQT0UtCr1BMRMuTiU8IT49NOM8HxGR0SkxpJXTZIxx7ynO5aKU6z80dX2UeHK04gpliXo/5XEdaOT0pKlIbPZkSQ92Kkffj0Xz2jtlHU8J7xMgho3QPiakn8gVdZlK4lvo8IiNCJxrNQZMTUx+pQukcEpNznyEnpoNuYZPvNXeqQj1kFJu2bvvl+WRqmaCjah3a8s7ZBPbYjPd0lvKd6MUmvPEFWUrk3718vvLiEAHmLiEiDYMbomWZEZPvaidkxVnnRJoTiYt2oZROuPWS+GVCwiXncXweEN5uVbguBESgFCCVv+9pTHg7tWn7h2pnbC29ri+vVlbXXxLCu1++vkhSwdhHpSCnVnAzK2Lnte8Nq/Al5aQyM1/a5XstwnHGJ3x6Ekrt+yTK4pPK4MYhXG/ZFyRay2nBn70U5PVMaH/tMj9NGss95xRWuRz76dK69DJZOT8RvbQY3BAtQ+KifoKHcIUXSxk9K6bOyBN03WBCC0jkCXpBoPRTQLguJeSIHlRYi/Opi6WUL2ZlNZZnvVB6b0lGBPqW+jzlwUeN4ObemHF3RvGK6qUXPZkP1tzf1f7bVAY3VjF0a+GdGd18oFSWsbbcIyN7rT6PNy6nERGVYZ0bokblI/Bf0GujuOHoXqwJcgtafyUHG6C8M1DZGeoeJybe1UIVXT6KwMd6Vl8FqrXYNWBt+/sx2GcMJs8GEF3jKivZe1GjHs4BM2S+10o7zbDsNwYzzxpZeS8O/2vtlk7mQ/3QghdNGOnHxQEiomVhcEPUoPytCMb1OrP7tRP8GqZX6dpTP2gp3IlivDjUBfOiueoU2A6pcjiA+H05uNZ+yiAjByuZYPq1MbTq3DPMIEtEq8TghqhB6QcxY+D1ttp3L9ZB/ue0HFqa2VLK96tFYPVaTq1Qy3b5je9GkSj291Pl5yTiel9anV70H2DuGSLaWAxuiBpUeCY7fczPYomW3hvjlVY5sPZMhxwYKt45imD4fBDpiuCpgPgXPvgVJyaueGBlbENEG4zBDVGDWl6Rz4VupJBeqlrJL/L/NZNG5mc5uCQrWtY6wGhR4Y2H4OlWkP26H13Wfgx+NorRz4YxYO+AK9WL2L0AXL9hZENEG4/BDVGD2vbM12mJxheLbtJIr1GdF+U1ixzSe0dePEle4WnOGFDssP7GGFyO2Wd5OVTHLgd8US3ICnvQZu1C114LrL89DPdX05gOemDbKecjItpgDG6IGqRYVdiLQ1mMfFErg7B0N4yxb+Xwau3tglsOhr+N1l+nJnU/WvzfdrYXasUNFBOUUmusTA41Q5hncQQvG++vr4DUlQF0Xe1A6Cs3HD12qN02WBjUENEmw+CGqFF7euE+LR9NXfPA89eUdrqv8jCMgQ9TaJfNslf9eMpkh+uirCh8zYexv9W5Y1SII/JVslihd/Qt2Yx8Tgt27JCf+ws/Ig+NwTnP0wi+50X69dKdqdrSV3rRftKP9B4FizRKf6EyuUaanRPRy4bBDVHDTLBfCMPbqQ9nEXy7He1vjiB4PYrozQiCn/RrJ/8EHF9NwLWv+Abgfho12zv9nJmbns/Lx0l1WE+Pw9etBydJjPS5EKwOTvQKvZ95MPzMhdDVIdi2y8llLJ0uGOFNGP22gzj1sR/hm2H4Px7AYesAMu+F4O0pNeHOIPdMDs7JIn49Ygz+8SB268s4q9e5qVGuRJGt0Vorly+/Z9RYHaJ0psZ9plfN6JCx2vjJXqPuz9l+dL0Trn1XiohePjKZHxE1KjctAu+r810syKIcnxDTsuuD2Pmy1xSLUA95xGRGe+HOmFC7LQvea8yjirE7xvsXeFq+TrNwnJkQIb1Pq0hAePu09741IRKLJU0WsyJxySHMVetVuj0iJLMjJz4t72PKLGylzyzNJscW9p1Vr8x1A5ERk9rn1vuZWjifIizdNb5zWQZivSh7VaHOZWw2zHzjqvwutbqdIKKX1jb9H+3gQETLVcgjn00jkc4BO9rRtV/BercVKjxOInorhtRj/XFMC9r2dsF2wApzjbs1NWmfOf1jAulcK9osVliU5XziPFJfD8N1IYeuP6jo+q0Zc+n6Clmk4jFEro4jIvPgWC8mkHi/+hHZGtJ7O78/jQza0LHfAhMbahGRxOCGiJZUuOvHwAfDmHx1FLGvnPXr3DxPw+/cjVPXtOG+ADLfOOXjMCKiF4d1bohoCXGM9pyC/2YHvB8vEtjoXjHj8FGjZygioo3C4IaIFvc4jUSxDXoX2vcUpzREsVp414aINgSDGyJqUCOdThQwfU/vg8uGwZ51rG9DRLQIBjdEtLhXbbD36AN+hG4skVfmYRjjf8zCdn4U7v1yGhHRC8bghoiWYEb/p2OwK1mM/qFWnh1D/r4f/TYP8p/GEPnItu4tx4iI6mFrKSJqTDaK0Q8G4LtaQMdpJxyHVJhNBWTvxxC7HkQEDngvDMO5d66BOBHRhmBwQ0TL8yyN5J005vMq74D5d8vItUNEtM4Y3BAREVFTYZ0bIiIiaioMboiIiKipMLghIiKipsLghoiIiJoKgxsiIiJqKuveWurZs2f4v//7PzlGREREzc5kMuGf/umf5NiLt+7BzY0bN+Dz+eQYERERNbv/9//+H44fPy7HXrx1D27+67/+C3/5y1/kGBERETW7f/u3f8N//Md/yLEXj0n8iIiIqKmwQjERERE1FQY3RERE1FQY3BAREVFTYXBDRERETYXBDRERETUVBjdERETUVBjcEBERUVNhcENERERNhcENERERNRUGN0RERNRUGNwQERFRU2FwQ0RERE2FwQ0RERE1FQY3RERE1FQY3BAREVFTYXBDRERETYXBDRERETUVBjdERETURID/H4AY5QzrHdR7AAAAAElFTkSuQmCC)
# compute the VCI for the last pre-fire to view the drought over the time period
last_pre_fire_NDVI = pre_fires_NDVI.pop(2)
last_pre_fire_vci = computeVCI(last_pre_fire_NDVI, pre_fires_NDVI)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[26], line 2
1 # compute the VCI for the last pre-fire to view the drought over the time period
----> 2 last_pre_fire_NDVI = pre_fires_NDVI.pop(2)
3 last_pre_fire_vci = computeVCI(last_pre_fire_NDVI, pre_fires_NDVI)
NameError: name 'pre_fires_NDVI' is not defined
# view the non-thresholded VCI
showImage(last_pre_fire_vci)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[27], line 2
1 # view the non-thresholded VCI
----> 2 showImage(last_pre_fire_vci)
NameError: name 'last_pre_fire_vci' is not defined
# map specific colors to values to show the severity of the droughts
def remapVCI(vci):
remapped = numpy.zeros(vci.shape)
for x in range(remapped.shape[0]):
for y in range(remapped.shape[1]):
if vci[x, y] < 0.35:
remapped[x, y] = 1
elif vci[x, y] <= 0.50:
remapped[x, y] = 2
else:
remapped[x, y] = 3
return remapped
# define the VCI mapped/thresholded values
def showVCI(vci_image):
cmap = matplotlib.colors.ListedColormap(["red", "yellow", "green"])
plt.imshow(remapVCI(vci_image), cmap=cmap)
plt.show()
# view the mapped VCI values
showVCI(last_pre_fire_vci)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[30], line 2
1 # view the mapped VCI values
----> 2 showVCI(last_pre_fire_vci)
NameError: name 'last_pre_fire_vci' is not defined
Dataset-specific Functions: NBR and dNBR Code#
# this code creates the NBR for each image then uses the NBR to create the dNBR. It can easily be updated for other burnt area indices
def computeFireMasks(pre_fire, post_fire):
rows, columns, channels = pre_fire.shape
nbrPost = numpy.zeros((rows, columns))
nbrPre = numpy.zeros((rows, columns))
dnbr = numpy.zeros((rows, columns))
for x in range(columns):
for y in range(rows):
nbrPost[y, x] = (post_fire[y, x, 7] - post_fire[y, x, 11]) / (
post_fire[y, x, 7] + post_fire[y, x, 11]
)
nbrPre[y, x] = (pre_fire[y, x, 7] - pre_fire[y, x, 11]) / (
pre_fire[y, x, 7] + pre_fire[y, x, 11]
)
dnbr[y, x] = nbrPre[y, x] - nbrPost[y, x]
return dnbr
# TA Code
# run computeFireMasks
dnbr = computeFireMasks(pre_fires_numpy[2], post_fires_numpy[0])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[32], line 3
1 # TA Code
2 # run computeFireMasks
----> 3 dnbr = computeFireMasks(pre_fires_numpy[2], post_fires_numpy[0])
NameError: name 'pre_fires_numpy' is not defined
# this code applies a threshold to the dNBR to show the level of burn intensity (unburned, low severity, moderate severity, or high severity)
def remapDNBR(dnbr):
remapped = numpy.zeros((dnbr.shape[0], dnbr.shape[1]))
for x in range(remapped.shape[0]):
for y in range(remapped.shape[1]):
if numpy.isnan(dnbr[x, y]):
remapped[x, y] = numpy.nan
elif dnbr[x, y] <= -0.251:
remapped[x, y] = 1
elif dnbr[x, y] <= -0.101:
remapped[x, y] = 2
elif dnbr[x, y] <= 0.099:
remapped[x, y] = 3
elif dnbr[x, y] <= 0.269:
remapped[x, y] = 4
elif dnbr[x, y] <= 0.439:
remapped[x, y] = 5
elif dnbr[x, y] <= 0.659:
remapped[x, y] = 6
elif dnbr[x, y] <= 1.3:
remapped[x, y] = 7
else:
remapped[x, y] = 8
return remapped
# TA Code
dnbr_remapped = remapDNBR(dnbr)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[34], line 2
1 # TA Code
----> 2 dnbr_remapped = remapDNBR(dnbr)
NameError: name 'dnbr' is not defined
dNBR Threshold#
image and more information: Normalized Burn Ratio (NBR) UN-SPIDER Knowledge Portal
![threshold.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAABNcAAAFQCAYAAACcZpjjAAAAAXNSR0IArs4c6QAAAARnQU1BAACxjwv8YQUAAAAJcEhZcwAADsQAAA7EAZUrDhsAAHl2SURBVHhe7d09s+TKutD5XINxjcEBcyDitr1U7uCePhGDSwwRrWqw9t7W8A0uXqm8yzfA23s7wKrqCGJwMU4fF7u07L4x4GJh3BsBseZ5UpmplEpSSVlSvf5/52h3rSq9vzzKTGWmXj6EAQAAAAAAADDZ/+b+BQAAAAAAADARhWsAAAAAAABAIgrXAAAAAAAAgEQUrgEAAAAAAACJKFwDAAAAAAAAElG4BgAAAAAAACSicA0AAAAAAABIROEaAAAAAAAAkIjCNQAAAAAAACARhWsAAAAAAABAopcP4T4/lJeXF/cJAAAAAAAAmG5MsdnDFq59//7dfQIAAAAAAACm+fz5s/s0jGahAAAAAAAAQCIK1wAAAAAAAIBEFK4BAAAAAAAAiShcAwAAAAAAABJRuAYAAAAAAAAkonANAAAAAAAASEThGgAAAAAAAJCIwjUAAAAA+PHDfN/+Yn754x/NH/3wy9b8/sP9DgBADwrXAAAAADy5H+b37S9m+/2HfIr8+G5+++UXCtgAAIMoXAMAAADw5LRQ7ZP5vPnV/PqnP5k/yfDrz5/Cb9//TOkaAKAfhWsAAAAAntxns/n1V7P5/Mn4IrVPf/gcPgMAMITCNQAAAABo+GG+//6bayL6yXz+A8VsAIB+Lx/CfX4o379/d58AAAAA4JTvZvvHrfw38umz+XmzMT9RtgYAT+nz58/u0zBqrgEAAABAF32hwZY3hgIAhlG49nS+m+0v/vXivPloEZII24ZXuLOPAeCkdtyk8jkm+L71584fzS+3ftO9Rhoh5fr68bv5JUyzNd+fJi3z2WzcywzsCw02P1d9rmkB2y91jbYfv//i9o0MWwJWEuI+znBXcf9e3dj96h7i7vMWrsmB09dt1wmH6gD+IjfO3x85BfH9z1EC6RbefCTHIRT2aXB8hDvr37j+OdQP8+OBTycAmEcrbrpPlR+N+4S9V/QmquJ7SjsheDyf5lClAb73Bu1T08t6JaUhlprvk/jxe5T4vod+sa6RRhi6vnrIisXT/I37dG0/vm8baffT+SvtN02maVxj/lp3owz49Pkn81NoDST7wU3TeNHBdzkHR+1UNBH3u+bnB+L+gLuL+/fqtu5X9xB3n7Jwzd6YJWD9JgGreUz0pHnsqt8/fKrghsQX6mUuWgDA3fgRPxSq6H28L5/Vfx85lTCs0gDbXySz1TnzU9PrsjUN8Yv5o2bW3HenLTXf5/Djz99lDzqfPhvyWA9KH4r/UhWwnLhcItX1vP1NpmlM5K/1VuGc1tZrF2jId3KKOZ/MX/rz69MfzOdwrt3Cw+oHQ9xPnO9zIO4/qTuIu89XuCY3ye3JG/PpgHevPv3hp+ikPC7pb1S3/OX3E/sJAIDr+F7neGc3lIkbRTNrZ82gx1LzvVvfze+/1SmVT5//UD/VxsP4YWuebU10qEfQmkFaI8n92eP7tlXTyRdohLTwb6Fg49PPkoZ2nzUN/Yc6QW1+/PY7BSAXQNzHLcd98tFLu/24+3SFa/UrtSufN7+GfhX+9OuvZvOzVjeU/z1q6uzTZ7P51fcl8StvPgIA3JHo/nxmk4BPP0f3f7kfbn5u3hBPZuI+b6LptW+mqLmCSl2/peb7iLSrC/dRzw2aBj0gfSiuNc/0s6ZhfR9op8h1EhfGfZLr6ld/TTWu9aj2w6efzK+/bsznVibgk007/2p+bSWaG02U5ExcsNznyRH3ifsR4v5Tu/W4+2SFa3VfCdann81PUemnRu7PP+nNl0InAABuj9ynw337h/lttn46q/t/M88d9/txmmbemxm1eWrBLzXfR9Do6oKmQQ/NFo7Zgi/3xQnNQpLP5qeosOJT61r/IZn1cCa5grRGQUdHgZvVaKLUOh8xI+L+XPN9BMT9J3fjcfd5X2igJEpNOxw9naJuf29VO5fxojeY9FUL1aroYZyj6r5jl+VoPxGNeek61PPf+om0arEfL+r001dj/SV+zPfjt+abosZu10Xe5HHDxwIA0E3vQb9E8VZi6Xbi4/i//ClqmtV4gj2zT38ZMuNjfQodMs3r9Hz9fSret+5e1dc07ijdoIenda/T/qoGbnRH43cNi91Tf5g/xxP0HS9J6+k9vvECK90uyaB3L27K+iXs91ES99EM19eQH99/t/2UhfkfnR+6r+PfZ0gnam2yrpo8g1oP0z9LZsx9rER9p6nJ+QGvOZ9GIR1qxP0kxP0uI+L+DNuoy6niXWvfyt+238fWpOPy0SNwv3I/DrntuPtkhWutm6mc5ls5mUa9icVeqHLQuzpF/f6bPTnqOKLVlxtH3Rz3t6fVGOsvm+NPWVaXvzG/23HqiX+E1zx1v/XjRztKdBq5XXE7+K4nfee6q2MBALAkcWs7Zm4EU0mU/SaJK02sum9O+2z+UOeyzO+zFB5UicP69vXJ/Fy/InC05hPUqCnTmQbn27hPtfeF3uc1rdPxVkUZtx5b0w2us/h4FvKHTeh37GMtROnKZAya9Z6q2+Y+is70hi7PpYeai9PzbnuckNd9Nbh+UQYpdb+fkrqPZru+utnzYPtbM7Ok+9F+71foRtKJdv+7j6JrGc3vWoVxEzS3t3WeYcbzkrjfPGWJ+6ozfsgI9SjTt9HeB2SdqnjX+l3nbQvqmus8Lh99Aver0XHhluPu09Vc+/xTq68GOZmqN7FUJ3M3PeGafbUdqy5Sf9432wPLb+3URaPa8aeovfj0ZR2RkzkuOB+jMzh1SN+uudzZsQAACImnQ6m2o8TesPhentqh7Q9NxMVPTP2N89Nn8/OvU7uHkMSdZNLCPNQszVVGzHdkwvL7Nkpot51INxztY0kYT06Ez31PbTXf6qrl0e5nd9gP8/vQPmqbY78fSd1H815fx2RbB+ahmU2/+OunE0Xr3OgyV22j5nzSC+keE3E/DXG/14i43zB1G2Xu9j5wcl2a6zw2Hz2E+5UzIi7cctx9vmahWr28q2q5HEg9+FqK2r64f/zevPA+/Rx1itroWDVKRLTbA7eqLPa9QjhpWT3qaX81mxOR/tNPVf8SjU5eP/0clv2nP22q6tiJ2zWXez0WAPDMjuKp9p8U4mnHPfkUuZfXFQy09rH7OIcf3yWRO+ItX/rkNcqk/aJPbt1PsoLm581P07dLTZ2vJOr170+ffz7uK6rVZ8+pBGicbvi5UYGjuY8b/VlFx/JPv8b3RllXfYGSHF81+z21kVnoqi3SXOfPG7eOdj03sn2tCb43O8BvbJemo2R8zUCFqWbc717qPpr9+ur0yXyOzo/ejuCvnE7scrogTc6lkcfoSHxOnDOfB0Tcn4C4P1Pc7zZ2G7VQ89R9oFb3/zc6H92L+1U13ci4cMNx9zn7XNOD2NcpqtZk+yUuYW217dZpJbL7ST99/qmRwKgTEc1Xxep86/Oqr7146rKO6dtw6mnHB5/TPpmf6jubrMRvpu5XtLn+878a+T6PBQA8u2atl1bn4hpfEzJacS2GUZmiCbQ5g+2DxP09hU+8zv1ipN75hj6pfpJ0TfPHZi2iYe10w0/tmv5BnPiWjFR0b2xmfuME7zXuqc20j54jYRU+6TnYfPtjuwP8TeOc/CSrLOP/GmVwZ9rvtfR9tMT11abbv4nOj89xH1gq1Ci5ZjoRt4S4fz7i/vmmbGNjfbVw7Og+cN4LMPpxv1Jz3a+u6XlfaCAn6kZL1rtKg+Xw/7b1AbtZGvrJ/I35/v17Y4hPkVhv1fgff25UN/0s41XSl9UkF9WSjwJbHcOGC7yxXUuswz0eCwB4dq0noUedi4tWYm+UOEHfyMCPownu+OmtTQ/Eq/EjvV+flM0ZY2i+vvPlZkfIcZON5n2tqeOePapj7+rJvZ+trkN9HOLMwjXuqe0Ha9r/i+6TqhuQ5vJGnKM9ztvvsdR9tND11dD1kPZz1AeWkJUPq3G1dGK302+TO+MhdOs6Ob2sZ0Hcn8vQfIn7p0zZxtb6yjl7PF77vhLFvbNwvwrGXEg3HHeft3DNc6XBzSqtolG7qWY759tum0NcLzO+yNpV4/1Z2LgQW4mTyKRlXVRrnb+7t/Y0tn35qv4cCwC4BxIfowA5R98kXqMWQ+NJbgJND2itdvenGnyS/nnjMmjNpiK29sP0jmlqk+er/Z5oh81V58vL34tamYDv25BQ13Wol99+iVTtUvdUbarTbMYjZN7aDUjz7Wgp5+iy+338PkpZ96XdRjpxyGyZsdlqrjya5c5L4r4i7j8a7lcT3HDcpXDN+9SsvpisUcrdDkRV4qJxI5hQ2nyks0T9MuIbm16wuknxdt1EVf8nORYAcNuatUKaTQLOFNdi0KfnZ886pQbLJ8kXtWo/zLQuY+ZrOyAO38n6f/7Z/CzTbXT4eZnmFZ+0SeDgjKt1v4V7qjaP+dOvkmlp74sf2qmyb6Uw/Ry9xn5vCPtowevrDFdNJ7ZrNXTsk+Z3/QUCU831ooT7t+B5Sdwn7j8o7ldpbinuUrh2UvMkOKpO3B5a7YSbzRE1cVElMLy6GaI6b1kX1aoJ9v3P22i7JLgu8jiSYwEAd6/riaMksFKTWHUmXpv8/439dHmtfqZkXXxHx+c5Nd/mfcx2KLzRjOdn81mHhaoGxQn15pNnuXe6PoLkthi58j1VFv75J9k3WiskzrX2tFI4/VR8if0+0z6a+fryjvNAzX2gK99Yn6ukE73mvpSFy9rGWk2T2uuO+RH3JyDuX0dzfbtrMLb7ZVsgdsg8uV/drycrXPtutrb6qnYS2Dxwtg1yXH0xVFdvPs368dv2uO2zzOv771vzyy9/PH5FcTtx0eh8s90M8cxlzS26iI9Llts1wb7XCZfFqvo/8bEAgLvVul9oXyKSUajiqYuj5wTRKLb/6O3nYwSJ6b9vf2m8gWtS7ZrPPy1Qi0FMmO+PH1FmwD7pjpvrzCVOqLuuNUIi+lebaI8Pd2WBe2ojUyPzOtpQTfP9Ys+t+jeZJl6R4PgctdO5P3X+dt3+KJlL901snv2euo8Wvr4smY9cG7o+9d/NfXF8rVwjnei1li1L/r21T+LrXF/OkLw6ckzqWTUznM+NuH8W4n63k3H/HB3nbPs+IHFv9LkymI9u434lf1TjjzkRZN71sm4r7j5X4ZorHbXtglud+TXbizdvtM2q7XphNae1nQ3qK5MbZ5TXPomik0EuyEZ5jjhvWedrXsQ/zG9yQdj980vztbnKVg92n2NnVfVvvIo6GtyF9kzHAgAeRfttVJoQq+Kpi6Pu+zTRU/44rp/w47df6njuYrqsSkQyEJNe/XZc2yB+/Xy6ofm2EpXaP0/YnrqWwXIkYe7SCfF+1P5hfm8ltGe/p47q0FjWQfsGitaxkXCPCnmO0jRRn0LH5+ky+z11Hy17fXnV+oT90bhWumuizZNO1ExnvA+a6fXv2+i3aKXay27uk2gOn36umxgmaJ53zQznsyPun4O432nhjuyPztn2faCxc4/PlSn56GPcr8bGhVuOu89VuNb7dpAmrdrYqNb66Sdzzmth2yeR12yG6Jy5rLN1FDL1a9f2UgtX9X+mYwEAj0Li6WC/ppLB/Xn8zedY+yn/uSQRu2l1cj1K6x7647cxCeoReucrGbBNnMC9hPZr/1s0o/vju/lNE9qSoQhJ4Nnvqc0Mw/S+W7TpTVxj6cR2NSy031P30eLXl0w/MHvtJ6i7POLC6cRGE6Mxx1PHOaPWmmicd3JCLrRl94m4fx7ifodz4/4JY9e371yZlI+egvtV7Jbj7pM1C/1s3wjz8+dWCa7SA2Pbi//J/NqVQvi8Mb/6DgY7pv3sppVz7JhWXW5/LydOV3mOdc6yJunYD3qxarBo/aBPGbsWd1QoFZWqj9O1Difc1bFI2D4AeED6JqxfJYEXVyC2cdR24PuT+cveYDnmwViVeGzPu/k081Q81t8/m5/t+uh90H0diafvnpekMxqJxlbfTj3Omq8mcLvu2zbx33pYGJPxOxc1gnZMPypL8+M383ucy5w1fdOqjX7Ur1aVnrFpPvdNRf62y+nYN3b9hqaJMlOp+71B5tveD4n7KP366iHj11P8paSf3Tq5b1S1rcPH6/x04gTth+h9x1P3i/Y79Kfu63y85vV9VsuNB0Xc73fWfIn7laO4L87YRmtgff25Uu1792VDdd8Zm4+ucb+y44+6X9123H35EO7zQ/n+/ehSwyK0un7dpltr/XUWTgIAgPv3fRs1v/MZCfen+vG72f7y22XSBY11+WR+lkQ+SZBb88DpRDnXf4maqX7ejC0gAO4McR+34kpxV1/+MAZvC8VZfvweV33+ZBZtEgoAAK5Kay8EHU+XtbnGqNoNc2jUBpirryPM6ZHTiT8aNXm6msACj4G4j1tx63GXwjWcof064gWr+gMAgKv7FGesfugLoqIOjGVovrXs88TOwadqNpv68T16CxpuwCOnE5vb9unn7pc3AI+AuI/bcPtxl8I1pPv+e/qrqwEAwN2xbykbdbPX/lOiPl8W0nhR0Y/vhkoMN+SR04k//hy9ZY+WG3hsxH3chDuIuxSuIVn7NbgkLAAAeHRVB8hVZ8TH932t4VB1Fj+2k+QzffrJ/BRyfTQRuiWPnE6MmyZ9+iznIElgPDTiPq7vHuIuLzQAAAAAAAAAWnihAQAAAAAAALAwCtcAAAAAAACARBSuAQAAAAAAAIkets+1l5cX9wkAAAAAAACYbkyx2cMWrgEAAAAAAABLo1koAAAAAAAAkIjCNQAAAAAAACARhWsAAAAAAABAIgrXAAAAAAAAgEQUrgEAAAAAAACJKFwDAAAAAAAAEr18CPd5lH/9b/6J+3R5f/1X/8V9AgAAAAAAAK6PmmsAAAAAAABAIgrXLqI029WLeXl5Meu9+won7dfVPnthp92JS5/nt3JdLbwe+3V1Hay2sqRnVpr9dmVWui/csNpeOTZwbAAAAACIpyxcKyWD5jNn/cPaUKRzf4aO7Wq1Nts9WWBUyu3arNdbc+unxH5fRaJ8szGZ/fSc9uuVWRdlVYiVZXZflPv3UKh1leOZb0xRrYj5RmgBeoWHZQPD0MOJcr+We3g8/sqsRlzvqdPdm3K/NetVK/0jf68HNvTcY1IpbWyuxifdDKSaM1YRb5c397YSwx/HUxauvZcPeJXDGjq2pWSACw0g1DJBuTVfi73Z7wuzvnbtpyGyntXq5SbP7TdPSo9V9SnffZiPw8EcPvRfV+B4teOZmS+5LV0zxS2fR8DdqmoGrySH0Ly9l5IZkevdZj7cVw2p092bvWyLbmchmbpWykb+tpmmhTbUZwaHMn8ATrmlWEW8PW3ubSWGP5rnbhaa74y+z6F72El2Fnfr6NgezM5mgkVZ2KCIifYP1AQuezWv7mOe+0+3p/wmN2/5Nys29xmP5jpnynfzbj/0FDJe8XhmX/KqgG+/54kfcIItHG/cm+th13Fta8aisMEjk2kP9fgHSaNVJesyzvHT9tTpet3q/U9jY5mZTNI8h4PbxsZ2Cln3oSTP1GNS7TutpaGZQdnDfjkAJps9VkWIt/ObfVuJ4Q+HPtfwJKogGIIMGeEnl5tduPHc6l1lb7bVHdy8vj75ne+9PJHAuuLxzL64BFBduw7ADELNXc08RA/IVCbX/Fshd3YlsXIbRYjU6e5RtjGHj4M5SOKmkUFqbKcmeeYLTpopszUddBlai3jTmXsDcMotxSri7WlLbCsx/OFQuIan8hoilz4pcB+BW3SqthZuhG8aOm/iB3h2vuauyQrTmfaXTIn/vtx/CwXwqdM9nFDwL97r/inP9SoZMs1YVjUr4twggCluKVYRb0+7+LYSw+8ShWujVW2stcNCWxhd7s12HXc8uLJ9/Yw78avqmOOnlfG32nFiPM2L/bu7nfQc6+qXWS9PO1Y8bmPuaVvzlM4dXdXUMI1rx+5+Xc6rBBf3sSFtO6p26/E0x0Ndpbd9fOJp3XctXR1d9h5/X526Z166vGp/9/xu9mbtlqHrHF4S4TdAm9W63+16dM9ETD3PzzVlef4YdHf2Oe14ts2z3fVNPAtNHmt+/X3V8/Z15K5V+1u/SedV0HHNrmQ9oknSz5mmo/lE52bz/O07nlOvtbTrP3t1R4gascBsQh+qcn31Jf/DAzMZt3oYkT5dl5RYlhZX70e22TVrbABIMmesOhfx9rRbOl7nIIYvi8K1FN8087U2RePClQxnId83c3bH3rVgo33RD02rmUkZv5CMsl7UctFm7sLVv20mtyPYBEnrqp0r+mW6r5T8YTOeX9tt4DUD6zLz+kNYRx2/6tyxexXdtsm6NRajy9BCgwWiUgiMsn7dBRbTt8MeA9duXSZy00TcfDrDmByfl5WfVsmyGttdFQ50dXTpj/9RnwR57vrnas+rEgpt5L/7rlcchgICV2NKbyKt9bd/u6GzyeLk8/xMMy7vrOM543r4czXLv3Qvy7OFpe3rSM9Zdx25b5oSziurukb8svx+0JeFNK6PlHOmS8d85As3n74C8h6jrrWUOCbCNUeNWGAekjFx16de631CwXa49lKn6zEplqXG1aW8V7FM6Xa4jwBuwcyx6izE29OucbyI4feIwrXJSlNooZOc4rZKpe3nJ2p3LRnaoTz0viimTWubhsn4xc69He9gDjLotIeimq4svvZk+lLWVYOVr4WSmSJMVy27q6RbA5jtGiorzE47YwzrqOPrGLIenQVyPuMfL0emKbSD8H2UCZ6Bq73ntzffuLcMRpK2w7a/txO5bWhPI1439rvjKsTV8dFCLLs8t4/j8cI6yfz1GNpzIOwntwVHL2hwhWKiq5laKGAUZUfpW+nvHq7AoOpkU5btN0j2z5v+7YawnZHJ5/mZZlveWcdzzu2ub+LDtFCrvUyJDWFl5feOhaadV7p7voZrRKex54BOF3e8KlLOmS5H87HXip+PLNN9e9qEa21SHGtrF9gBiOl1prUK/LC2tUK7rqooUzFEMltR6BGp03WbEstS4+piopq0uU8UdBh/TADMZ95Y1YV4O2e8Xf54HSGG36XnLlyTzG5cjbQe+mp7eFUGry5o0gv7zfhr+X2wqHritFnV0eFOcoHtizXb+DcIDmXoJi5vvw3BqpDpNo0ccxXAPg5RwZQtjNAPMv6brE9jJXV8lwGWAFeN55TfQuGZzrNejkyz0eBad+KYpH1sfe29LJft6nh7SuJ2hFpgubazjyfS7XDbMNhUTPdzVDARP7EJ61TtIz2G9RJ0/nIz8QdSjltcwBqC8NGy5e/4i6PfS/PNHZihQH7aOddIinmWd/7xnGu765v46Vpe7WXKUjWBEs6N1vqecV711qZznaJuGl/emjHX2sQ4FujT1urT/Oc28EgkvSJxxA97VytU79F9cXV0TdeW1OmSnBFXlxE9WJEMatfDoNr0YwJgPsvFqunXNvH2tMtsKzH8XlFzLYHNaB1dV1pdtfrUVRvIO2faY3WGrs/U5YVaS5JZ/jIidjQKIzrHr2tRxZnOuj+p051CzkczzrvO9UzdjnPZgsWe/XxyH4nsiy90lYAaH0pZ2WqSVrVkGcmOpp1Z2i9ahW2h0LPe3hTznuenXXp5fWZbD3+cZNqB2udW9zJlynBuNI/xOeeV70uiLIb7ILtFo6615Ov/8ucacE/yndZ09TUL3KAPBXwGyDYtv2STyXmddb9ewH7tM1XVA4OusPboxwR4VsTby8bbJRDD79dzF67lu+YJGQbJrLpRpgidGCY4PW1p9tutWTeqfromWgn6ljems8ZYGH+/jtarOfiC9zjTOXU5k8XH9rBzNYdKU8j6dD3BSN2O0Lb+6MlIabZbba4nQkHXNKP2UfQmmWam3xcGyHkT9avmb0hZvukuLJDP9i85P06U6SQ55xpJMXV5Sx3PS2+3lb12HsNzzitbW9Z+rzds7TRWbv5+fpPtqz7Ouoblq3U0pF7/wHObdg0fhUH5QmsYfPiq5D01Q1MfaKVOl+K8+3WX9Pio/W3WFR7eeh8qqNRjAmDItOt3iVhFvF0u3i69rcTw+0bNtTtQvT1FMndFYfb7UjJ39TCvcZ019onXKx68lHnOIsvN5nCoC9hO9JvUtQ06eI3tyDdRwZ02Qa2Csf5ru3iSsF8kVcEb2+dWP1+gE78OurohSeD9IoMrXYt/9320nexE/1EtdjxvxbnnVbMJqn+Rgb5JuOqrbpqu66warleA1b0+Pdc/8OS6rpVqmHANRw8s6kyLZJrG3IRk/GoK/4KT1OnOcf79ukv3ftVhYGH6gMDlyrLiYA5DubIhnccEwFjd164O/vq9RqwSxNtO3cdKhwsfL2L43aNw7dbpReaqp1UdPcY17HyB0VzOaNrUWwuwHpIDxCwys/GFIuXedL0k05q6HXvXR50GMRt1q2Csy8ty7Rg9tR+q+likClWiw/b65oEu4PvAe/R7Vfj2lBY7nrfi/PPKnh+7g1wLctPfyTlmd1NpinVdu2uc3Ow6ri87+Cdul3bzcQy4JZe5hofSI6Eriw6p0003R1xtS9i35bbOlGl6kVgFXMm06/dyseo04m1ruOTxIoY/BArXblx426NeZHKBjyo1n4MEB1veckJo7jZyfC91urOEUvzSFK06sqnrUx0frc20swUvdUDWgod2x+iJhtYpejHEUQeboUp0KcFe/pF1tVsd9oP8az+43+U/9pYwsr+9R3SR4zlGaNLpjk0KfzxlezorW6WeV4EWOGqBv9Zkq77pejvtPTg/Hp1X6xdAW5SBGrguQ/OgcF9LnW4mZ8fVRJops2+OFi69COCWXTlWNRBvT1t4W4nhD4PCtTvRnWEbCCqJwhsi+2p32Yu/blYZ+qgaqg3WoZ6upy34fr1AG/HcbHxVP1/Q5KRth89Qa79mrnPNGdXHor+9fOjUU7bNj17LzBdXGqSFHv5pSvwWUP/Z/h76YzvRJFRuHKnlPbdt2eOZ6lR17v22+21Ae99HXKtq+vnnVVvdxLg3wXHj50xqHKvJdrvpLvMWKeBxlNuti2HNWtMnY5WkR/z3Xfe1qdON0hPLTi5TTIurU2jzfBfvs2KWTFnfMQEwn0VjVQ/ibbrltpUY/kgoXLtxod+s9tv5JBO4tp2Ju7/nIsGguqS1z6lmP0rawWIoVffiPqq+dq9Pud9W6+r+tsJ0ElLW8UsG9MUNK/OyXqZwI36D4jbupDJpO+qnGGUh++ZF++hyg+w72xHmWo9b4pa09tG6eQLY70KT4WLjjltT2N69bK+dvnWTec1avw8EYT+u7AHt+8+SbUvdvDY9v6r91zpXLmbh4zlJ3bfDUPVzS2PBS3ytVteQb6J5dG4kn1el7fB1LSmH5h7Qlz24hckObJw9C58zs0mNY96pWoLAE9uvq5efaIxqXFoSCOJ4o9dhoxVMK1Y10yPRU/72G+NSpxtyKpYlx9VzaVx2cUkzZYfut8q1JR8TAPNJjlV6nfr0abMfaeLtgvF2geOlvxHDH8zHRH/11//n1Ya57HLzoZs+PGQfxcFNYB0+5IKyv+U791XLociqaY9GOGfa3Ydcn9VvOmTZh1wj7m/5LH8fz/ec5YlDES2jPWQd85R1dMvz4+h6ZY3vchmrZZf3LyeT8eV3+7lvIzqEY9s7Tb1vJIjJX7GU7YjmF8Z3+zYassbJdPr4BIPHohqywZnE6yfD0Taf+j3WGtcN9bad3q6h8y781nWudDpneX3Txts43/EcvN56hGk6j0m0zHzoOuo5nknnVbxvZDh5bajWNG5o7r8RfCzoPTf6jsHpY1NLjGPq5PoBz2tMmkvjTVKs0rSCG7UhdbpeI2LZ2ffr6er75oghWvY5x2T0MmfeVuAhJcWqOG/Y/J14Ww9zx1tr5uNFDH88T1lzLfSvcxdyszvsTOE7e9KSe/nHd67+NltxfCTbmINbZryn/DLlgm7RdZTvC98nnD5dqJ4wSObU5IX2y7Q7fnJg+2uS7xvHY2D8WQy92CBhO/a+A3ztCP0gv1eDXFv2753cHZSteWg/TaTHws6n3d+eHJusWt/h6sN101B13OTz1O8x2Xdvcl7EKyLrMHt1Y5n/1Sr/LH08J6ibKp5oUtl1HYXztecJWNJ5Jcdfvg/TuGvDTqPr8NF1zV7onJlFYhwToQNbGY+Ka0CTvgDlsCtsjGpe+Ro7qnSFxpvOqOBjVSs9otfaYFohdbpeI2LZ2ffr6cJ9YqKzjgmA+STFqte6FYqMF0cB4q2sw0Lx1pr5eBHDH8+LlrC5z6P863/zT9yny/vrv/ov7hNwC7Q5nlbllSDeWbCg/DiZKSTQUTW3T9XkUAu2snNePX2WWzuefll6E/2QG3n1baXeX8e/4bJu4dwFAAAAcE30uQakar0UoVPoiwmDfI2xrDBv1yqcuLnjWfePd69v4XwK4Y1U2Y3WyAMAAACwNArXgFRRh5vaobvNX8dK+f6r79wyN+S7+5SuQ/zMFG/jOvJcxA0ez3xTVOs0puAPVxHeSEWnsQAAAMDTonANSJVJZjqvctPh7ZL6Rkkd7Oe1q42lbf+vWGh061yttax4u27hxC0ez+yLqVZJ35hkv8FNKc0397aoya+YBwAAAPAwKFwDznDUoaTtAF0G+agdahb6++SOPJ+MfXnAx030VXV7xzMLL+DYb9uv78bVRc2ZR79iHgAAAMDD4YUGAAAAAAAAQCJqrgEAAAAAAACJKFwDAAAAAAAAEk1uFgoAAAAAAACg8rCFay8vL+4TAAAAAAAAMN2YYjNqrgEAAAAAAACJ6HMNAAAAAAAASEThGgAAAAAAAJCIwjUAAAAAAAAgEYVrAAAAAAAAQCIK1wAAAAAAAIBEk98W+h/+/M/dp8v7F3/4j+4TAAAAAAAAcH3UXAMAAAAAAAASUbgGAAAAAAAAJKJw7SJKs129mJeXF7Peu69w0n5d7bMXdtqduPR5/iDX1X5dneerrWzRMyvNfrsyK90Xblhtr3xgOTYAAAAARnjKwrVSMnA+89Y/rA1FOvdn6NiuVmuz3ZNFxm3Z76tIk282JrOfntN+vTLroqwKsbLM7oty/x4Ktcrt2qzXW3PRSzjfmKJaEfON0AGcpdxvzXrVukfL3+uBizo8ZBsYhh6ulPu13Pvj8Vdmdek4clGljaXVto5Lx563j6YvD8CxuWMV8XZZy20rMfzePWXh2nv5gFc5rKFjW0oGudAAQi0U3Ipya6rKWbnJc/vNk9obV8Zo8t2H+TgczOFD/3UFjrKfvhY6TmHWF63NlpkvuS1dM8W1a9EBd2svmbwXSewXkthv3X3lb5uwn736cVWzeSXzbS6ylIyIxBGbyXRfPQifmR7KPDedt4+mLw/AsbljFfF2WcttKzH8MTx3s9B8Z/Rlqd3DTrK7uFtHx/ZgdjaTLMrCBiJMtKeJ3NzKb3JDlH+zYnOf8Wauc6J8N+/2Q08hY/ZqXt3HPPefLiP7klcFfPs9T/SAFHp9l5nJ5L58OET35YOks9xtWWPJ0G3ZFro37un1sOuIGZqBLGxQymTaQz1+WKZmMic8pb/p+1+VYa4y07LFfp+ekL6P0pYH4NjssYp4u6jZt9VKi6np65K2PIxDn2t4ElXgCTcFMsq4ur3ZVndF8/r65He29/JEAio3O5doCIXkl5J9cYmUunYdgAmyjTl8HMxBbsCNRHwm1/VbIXfnim8if7ZQI1gzidGDNdVYpsTg7XDkuQeaSbI1D3TbtNbvZsSjmjP2UdLyABxbIlYRb5ez0LYSwx8LhWt4Kq/hTqNPdtxH4BpO1dbCjfBNQyVBQukaMK9QeC3e6z4Wz+FrBJusMJ15Bsl8+u/L/bdZlnlNr5JB0sxVVVshymANOGcfpSwPwLGLxyri7VmW2lZi+GOhcG20ql2zdhJoC4DLvdmufQeA1ffaF9C4C6mqjjl+Whl/q50VxtO82L+720nPsa5+mfXytCPM43bdnrbvTulQ0VVNDdO4tuPu1+W8SnBxHxvStqNqtx5PczzUVbDbxyee1n3X4tvFx/PrPf6+OnXPvHR51f7u+d3szdotQ9c5vCTCb4A2q3W/2/XonomYep7P75r77ZT6xpiFJo81f4746tzt68Rdi/a3fpO2P+i4JleyHtEk6edE09F8on3Y3M9+f7T37dRrKe36zl7dEaLGK3DzQt+rct32ZRvCgzYZt3rI0S0l1qXF3XTZZtestTDCOfsoZXkAjs0Zq67lmeLtUseLGP5YKFxL8U0zZ2tTNC5cyZAW8n0z53fsXTPo7Yt+aFrNbMr4hRZqyTRyoWTuYtG/bSa4I9gESeuqnWH6ZbqvlPxhM6Zf223gNYPrMvv6Q1hHHb/qULF7Fd22ybo1FqPL0EKFBe4iIRjJ+nUXaEzfDnsMXLt1mchNE3Hz6QxjcnxeVn5aJctqbHdVeFDNv7lwf/yP+iTIc9d/V3telVCoI//dd70CMRQgyHx0Rhq4W+tv/3ZDZ5PGyef53G5gv53gz8Us/9J7Y7RsoV77OtFz0l0n7pumhO23qmvAL8sfd30ZSOP8TzknunTMR75w8+krAO8x6lpKiVMinBvUeAXm9V5dj2ogkzCeZCLcda8xpE8oMD91TU+Kdalx99Jm3kcAElzjOiTepruluEkMv2UUrk1WmkILneTSt1UqbT9AUVtnyfAOlR3si2LatLbpmIxf7Nzb8w7mIINOeyiq6cria0+mMGVdNVj5WiqZKcJ01bK7Sro1gNmuo7LC7LTzzLCOOr6OIevRWSDnCwbi5cg0hXYgvo8yyTNwtff89uYb9xbCSNJ2lNrm3U7ktqE9jXjd2O+Oq+1Wx0cLY+zy3D6OxwvrJPPXY2jPgbCf3BYcvaChLtzpasYWChhF2VGKVPqI7QoUqk5RZdl+g2T/vOnfbgjbGZl8ns/sFvbbsPrGOEwLtdr7Uq79sNPl946dmbb98tX2a7gGdBp7jHW60DlqJeWc6HI0H3st+PnIMt23p024libFqbZ2gR2As0S1QXMfgDvo9au1EfywtrVNu67WKPM4RDJpUUjrNSXWpcbdy5t3HwFIcYXrkHh7hluKm8TwW/bchWuSyY+rkdZDX20Qr8oA1gVNemG/GX8tvw8WD0+cNqs6ptxJLrF9gWQb/4bBoQzfxOXttyFYFTLdppGjrgLYxyEqmLKFS/pBxn+T9WmspI7vMsgS4KrxnPJbKDzTedbLkWk2GlzrTjeTtI+tr72X5bJdHW+7SdyOUJsp17bt8US6HW4bopvZMd3PUcFF/MQmrFO1j/QY1kvQ+cvNxB9IOW5xAWu4aR4tW/6Ovzj6vTTf3IEZuvGeds41cqa72G/1jfF0La/2vpSt0ARI2IbWupyx/b216eS60XXYNL68NWOupYlxKtCnqdWnRc9d4KlEDwckE3X8ACom6RyJT37Yu9qmem/vvGTF6Bq0czgj7l7TRfcRgE6XuQ6Jt3O5pbhJDL891FxLYDNiR+eyVletPnXVavHOmfZYneHrM3V5ofaNZKa/jLheG4VLnePXtYHiTGnd39TpjhjnoxnrXed6pm7HufSm0FdgcXIfieyLL3SVG2B8KGVlq0laVYFlJDuadmZpv2gVGoVCz3p7U8x7nk9zF/vNz0/WYqBGt9W9L+NtaK7LOdvv+2coi1N9Jd6eUddS8vV9mXMXeCb7dVR7/e24NrnKd1qD1tdIcIM+bPAZJ9tk/fpNLc+67wDAwoi3wGU8d+FavmsGkDBIZtaNMkXoODDB6WlLs99uzbpRVdc14UrQt7zQ9E1ykmO2Joy/X0fr1Rz8g5I4Uzp1OZPFx/awczWmSlPI+nQ9wUjdjtCe/ejJSGm2W23OJ0KBzTSj9lH05p9moYAvLJDzJuofzN+QsnzTXZggn+1fcn6cKPNJcs41MtYj7rdO2Wvnss7Zflsb1n6vCagXOe8lMebnN9m+6uOsa7hwtY3U6xvAkPRrXPuMrCtRvPUWjKuj24Z8oTUTPnwV9J4ap83Yvqzz7jtdLhM/L7mPgOcx7fpd+jok3p5at9s6XlMQw28PNdfuQPX2FMn8FYXZ7yUTL0HFD/Ma10Fin3i94sFLmecsstxsDoe6gO1Ev0pd26CD19iOfBMV3GkT1CoY67+2CygJ+0XfY5VBY/vk6hdqIUWvYK5uSHKj/CKDKyWKf/d9jZ3sZP9mPft+O3f7m01Q/YsM9E3BVd+C03RdR9VwvQKs7vWpt+1qcQq4Q13XUjUMXONayO1yellxMIehnN4QicX+7lpnMCSzNWZ2Mn41xcQXp3Q6/77TpXu/6nDuwq6xj4Dn0n3t6uCv3wtdh8TbUbqPlQ4XPl6jEMNvGYVrt06DoqueVnX0GNew8wVGczmj6VNvLcB6SA7os8jMxhdylXvT9bJHa+p27F0fdXrTsZGuCsa6vCzXjtNT+6mqj0WqUCU6bK9vPuiCrL9RHv1eFSLdp2ffb+dvv92OXdXRv3YqW53WpSnWde2ucXKz67h+7HDU6eGF3HycAu5JwjVebuuMnqZpFrzehtIxoQuMWcwRd9suEz8vt4+AZzLt+l3sOiTejnQjxysBMfz2ULh248JbCzUoygU+qqR6DnJB2nKDE0Izv5Hje6nTnSU8dSlN0arTnLo+1fHR2mk7W5BWB2QtmGh3nJ5oaJ2iF0McdWoZqkSXEmDlH1lXu9VhP8i/9oP7Xf5jw/DI/vZu3i3vt9Ck080jhV+unH+dla1Stz/QAmIt0NeabNU3XW9RvQfnx5v6SSk124BEmtGzbz8WLk0zvyjjNXC9h2ZFIa7P5Oy4ewlX3kcAxMLXIfF25nh7S3GTGH7LKFy7E90ZuoGgkii86bCvdpcN1nWzytDnWN/4Perpet7Ot193f3+W3Gx8VT9fYOKkbYfPcGv/XK5zzRnVx6JnH4nQqadsmx+9lpkvrnRPC0X8E4z4bZb+s/099Ct2ommjBOvU8qBLuNn91uNUfwn7bffbmfa+T79Wde/zt7+tbgrbexO/8XMiNU7VZLvddNfNFAP3SpuYu5iVFbNk9Mrt1sXGZq3hkzFQ0jH++ziuj9IT6+aPu8tadB8BGGW565B4u0S8XXRbJ7qldUEThWs3LvT/1H57n2QS17azcff3XOQCrC5B1/l/tADtEDM8BfHiPse+dq9Pud9W6+r+tsJ0cgtYxy8Z0Bc3rMzL2gfFedVvj9mbbdxJZdJ21E8OykL2zYv2ueYG7QRTh7Uet8Qtae2jdfMEsN+FJsPFxh23prC9e9leO33rJvPqXl8dfh9o2ujHlT2gff9Zsm2pm9em51e1/1rnylS3tt861f0lDFXptvRaf4mvxeoa8U00j7YheftL26HrWu7G8RT2e3+HlhO+sZULnxOzSbq+I6dqCQIYoLHFXVua0Tt0v6mubb+uXqqisa9xyUqAieOYXt+N1k6tGFjHzuo6D+mYgTfNHTkV62a471zUEvsIwDTJ16HGFJ/faPcjTbyN13fWeLvI8Uq0xH7HPD4m+vff/++rDXPZ5eZDN314yD6Kg5vAOnzISWx/y3fuq5ZDkVXTHo1wzrS7D7kmqt90yLIPuZbc3/JZ/j6e7znLE4ciWkZ7yDrmKevolufH0fXKGt/lMlbLLu9fTibjy+/2c99GdAjHtneaet9IwJG/YinbEc0vjO/2bTRkjZPp9PEJBo9FNWSDM4nXT4ajbT71e6w1rhvqbTu9XUPnXfit61zpNLC8m9pv3cL2dk4bbVs+dJ30LDdp+9vbNOIabk/jhub5PoK/1nuPfbQ/GiP0fd8lMU6pk+sHoE8d20cM0YU8Jq2mcSwpBmoaw406zohYd/Z9J83o/dteduI+Sl4egGNJ12GcN2z+TrythyXi7dzHSxHDH8tT1lwL/e/chdzsDjtTuGZqtuRe/vGd5b8tURqdbczBLTPeU36Zx7WLdR3l+8L3CadPF6onDJJ5NXmh/Tbtjp8c2P6c5PvG8RgYfxZDLzZI2I69f6GBdpR+kN+rQa4t+/dO7g7K1jy0nybSY2Hn49fJk2OTVes7XN27buKojpsunvo9JvvuTc6LeEVkHWbvxF/mf3bloJvab93qpoonmlR2XSfyuTofe55IJm2/HF/5Pkzjzn07ja7DR9c1eaFzYhaJcUqETmFlPCquAdOEWDeRvljlsCts7GtGFI1JVXpE41hntPExsJWO0Ws4LY0xItadfd+5sNn3EYDJkq7D17o1hYwXR1jirazDkvF25uN1FmL4TXrREjb3eZT/8Od/7j5d3r/4w390n4BboM31tOq1BPHOggflx8lMIYGeFxH2qZokakFlds6rwu+KPzc0UfMhN8fq20q9P45/w2U947kJAAAAYAr6XANS7ZsvRegU+mrCIF8DMCvM29MUXtT9uN3rWzifQnjjVHajNfIAAAAAXBuFa0CqqMNN7fDd5r9jpXz/1XcomRvy5X1K12F+Zoq3cR2vPop8U1TbO6agFlcR3jjV7sQXAAAAABwK14BUmWS28yq3Hd4Wqm8I1cF+XrvaWNr2/7kKjSZxtday4u35Ci+yL6Y6hfSNSPYb3JTSfHNvYOJV5gAAAAD6ULgGnOGoA1DbQboM8lE71Cz0dzqUHGZfBvHxpH1ZZeEFG/vtTK/nxnyi5sq8yhwAAABAH15oAAAAAAAAACSi5hoAAAAAAACQaHLNNQAAAAAAAAAVaq4BAAAAAAAAiR625trLy4v7BAAAAAAAAEw3ptiMZqEAAAAAAABAIpqFAgAAAAAAAIkoXAMAAAAAAAASUbgGAAAAAAAAJKJwDQAAAAAAAEhE4RoAAAAAAACQaPLbQv/+f/2n7tPl/Y9//J/dJwAAAAAAAOD6qLkGAAAAAAAAJKJwDQAAAAAAAEhE4dpVlGa7ejEvLy9mvXdf4aT9utpnL+w0AAAAAABwI56ycG1KIQ0FOgAAAEspJa21qtJaL2szJrVV7tdm5R5SVsPKrNZbsy/dCD1Sp3sIZbWfV2Hb/fbv5QgMSJ0OQLLlYhXx9lJm335i+F2g5hoAAAAurtxvzXq1MuvRuY2q5r/NFDQmKWVehZuX+6ohdboHUW4lk1ftZ938LMtMVv3gMoA9mezU6QAkWi5WEW8vZYHtJ4bfDQrXAAAAcEHuSfq6sE/xJb0/ik5T2MxKZvLdwegL7+1w2JnczkPne5xZSJ2ul2RK7NP/1dZmWG7bXjJzRbWe+c4cZLsPh4P99+NQVBmtUsY5yu2lTgcg1eyxyiLeXtL8x5AYfk8oXAMAAMDFaObD1p7IcrPTxP4md78MKLdm6/IAmmHZVbmUis7nzWUWJEOx3UZZsNTpHkS53bpMnGzrLnfb6mQbc5DvrL3sp8ZuS5sOQKKFYhXx9oIW2H5i+H2hcA0AAAAX8yqZDM14VE/yG0n+XuW3ffUEPitMZ95QMgv++3L/LdRwSJ3u9u3N2vadM1QDojTfNFOtctnn1aem8H1pynf7QaROByDVUrGKeDuHMfF2ie0nht8bCteSVe2ptUNAW9pb7s127TuIrL5fb90FdpKrrhtNa9tpu1+b/MXd/6bR/pcwtNdZ2967cf131hzb5tty+2mq6U535NjeF67NuvsVAADct2yzaz7RH+Hdd17z+tp8Ah959RlHGdfnFVKn61JuXVrIp6/Kople6Xj87/s58uPY8VwfOJeUZa/uU5vsF7f57+/H65Q6HYBp5oxVMeLt5eLtUsdQEcPvA4Vrc/hWdQhYNC7c0uwL+b6vBMx730qgqC76empfMDWlPfZEss4vq6rtfaWn1Dpp27RwzhUQ6mRy5WoHitV2VR05dsRDoQWH7X0hU+m+0CcFUyIQAAB4EJIJcWmA/oyC/CYZmsq7ZHL039TpemiGSdM07k9l/3bD62v8i+/UWtNazZmW8rc+SLyNPoR0vatP5aTqC6nTATg2c6w6C/E2zbWOITH8llC4drbSFIXWrIo7LYzaWO/XvTXM1L7QjgabHR6GNtCLdTJYrbMxudkdqmV+HA4d1VfTtk0DmO3IMSuq+Wv7fhmqaXUMme/XdoDToOgLEzNTxMsrtJ34/kSNNwAA8JjeJfHvPg6RzFac3UqfrlumnUJrmsan0ySd8+bSODr4r1VIC7k0lO1EOqRr3NK0JsYi6TwvyjyNaYIkOcNqnNTpAKSZN1adh3ibZoljSAy/NxSuzUILqfQi95eKXthvxl/Lw1Ut29PK1BpM/MT7/UK11zJThLeW6J99F/rEbQsdOcr83zb1/C2dVpapHyXA+Q4frfJbKDzToLiJl7fR4Oo7gAQAAM+qWVthvNTpknR0al0vXdM1kvkL6bxlO5POc5cDlXTX185mVOtmesxJnQ7AeS4aq04g3qaZc/uJ4feFwrUZ2AKjo2toXFXL7mll6i/+rR57LV+bnS3AGnHdT9220JFjvumZf679J1pxwdyUDiABAABu1ck0jajTeT3dcjTU/e3Wg6/tP/SbyLW/pepjWVR9EWkfRL5PXNuFR/Wz5gjdOonU6QDggm4q3i6BGH5XKFxbUOiwMEX2alz51U3q27bQkeNe+2rTC/h48DVy44K5MR1AAgCAW7Wv+lvtGmasKpDa8fIlO2welabJvoSHl0uvW76ruhzxSTftg0hXMcty2w2Hz4C1pU4HIDYtNl4yVp1CvE0z9zKI4feDwjUsprqAjwdvqLNHAABwX7ru+dVwsqrACZJp6s01RSRDU6UyXrVbG/tv2nTnqDu1nk9udqEPIT+4bjYGf6tV/Rc1xzscdlE3HN3pstTpANS646IOPlhcI1b1Id6mxdtlt58Yfh+eu3DtZOd94976gQ5yIccXcddwGNMuFQAA3IGuTIcbZnw8PlRQVw7kslKnm67uOuM+1N2PTOsnKHU64NlMi42Xi1WnEW/TXPYYEsNvyVMWroUmjWVphk/v+q0fFz/p5KKs1i3TFqJ3I+xbCRzDBZdNqdMBAIBHF2WgBtIJoXlQnruaBKnTzWQoTRO9yOmqGRvJXVX5q7pP3FFSpwPQ4cqxqoF4m+ZK208MvylPWbgWvyxg7TsA61But4ufdPttdyeI+23hLsp2ldG6yuneFzfHojemXEPmo0q5N996I9yxerrWW0Q93mgCAMDTit981pkeiN8YFyXaUqcbRTJJXQ9pTy5ThE64r5mxkbSaTwdnxWZ8Ri91OgC9Fo1VE6Wuy6LbcAfxdtHt70IMvznP2Sw025i38EretVmtt6E025KLd79emVVRfbnoSacn98tKLjS/ArLsbd3p//GyM/PFt5GWdV83ptOXCPhCuSvJN6bataUpvq6b+9Up91uzXrUKFcN0ulmyP8J01f54kR1y1e0CAADX004nRAkMTVeE9E/7jXGp0w15zcJD2r2fn6QdfYWE9jLrtJq6YBpTSWZuLbk57d+pJutqXzzl0mKy7W/trjpSpwOQJjlWaUxxb65cbatxzkW8TbPEMSSG35ePif73/+//utowr8PHLs8+dBcMDVlxcOO3HT7k4rHj5Dv3VcuhcPM/GsFPm33kRf4hp3RYXmPIChmzy+5DrseeafKPwm9X73L717lyzrYpWT83fTVkH1mmQ/xdLmO17Ib2hYwvv9vPwysPAABuWEhDnBra9/tD0Z9O0EHTCm7UhtTpetXppHhopBlPLVOG7BLpGZ926hv60pqp0wFIlxSr4nzh8e/E22q4SLxVcx9DYvhdeeIXGmQSQw72bRl55kvEPfk7L8zusHyn+9kXfYNHtQ6BfM4L/X7TWi8vl3XbmSJ6y4fdHjvNznxZdpVH0PU7mF3hX/2rTxeqJwxZ2LbjN1nJATneF9F2HY0PAACeR7Yxhw9JX0j6p5HUGUpbqNTpemVm8ybpsDi9kuUmjxNgfpkhLeTJOsi4mk46zPiih175pmMddNNzU0g6+KMvrZk6HYB0SbHqtW7qKOP5br/ORrxNM/cxJIbflRctYXOfR/n7//Wfuk+X9z/+8X92nwAAAAAAAIDre+KaawAAAAAAAMB5KFwDAAAAAAAAElG4BgAAAAAAACSicA0AAAAAAABIROEaAAAAAAAAkIjCNQAAAAAAACDRy4dwnwEAAAAAAABMQM01AAAAAAAAINHD1lx7eXlxnwAAAAAAAIDpxhSb0SwUAAAAAAAASESzUAAAAAAAACARhWsAAAAAAABAIgrXAAAAAAAAgEQUrgEAAAAAAACJJr/Q4O/+1z90ny7vL/7ef3efAAAAAAAAgOuj5hoAAAAAAACQiMI1AAAAAAAAIBGFawAAAAAAAEAiCtcAAAAAAACARE9ZuLZfv5iXFxnWe/fNgyj3ZrtemZVumx9WK7Pe7k1ZunEeVLldm/V6a/YPvp0AANy3UtJhK5dOWZsxKbFyvzarVZS2eZG0zoh7fup0j6Dcb81a0oD1tsugacKRG19KmrJr+pWknTvnIAlNPa6NNKjd3z3jAwiWi1XE20uZe/uJ4XdK3xY6xd/+z39wtWEuu9zoG1I/TL5z39y/wy6vtqlvyIqPgxv34RyKj8xv5wMdUwAAHslhV3zkWZQ2MfnH8F378FE0xm8P2Uf3bT91ukewa+3jjmFw42Xf5Vlj/CzLZPB/Zx9FO0EZp8P8+NHfJjt1nIFntVysIt5eytzbTwy/ZxSuPYL4gsiKj110wRwOu+oCe+TCNQ1CbvvzeOOVL3R86O0HAOCWHeR2XCf260T+cII9pNds5iRO3MSZj+N5pE7X657SEjZNKBkjSeMe4pVtbLuml9z3LfG+K9ppKnHYyXzd50qdBtN0deM3Hn4Cg2aPVRbx9pJm335i+F2jcO0BHAofQKcG3ydA4RoAAFcV0l369FtzC6G2/UC6JUrUdybXot+z+DF86nRDHiUtcSqjFI5LR82GHifToAnzBJ7CErFKEG8vaIntHxLNjxh+m3ihwQN49x2qZZl5rT4BAAAsaG/Wtk+W0/34vGa55AMO5uOwM7mkVcYov7l+XrLCbHL7VVO2Cd+X+2+hT5jU6W7f+P3dK/si+999fn9vbXtptttqzlnxZjajDlNpvvn+f3I5xtWnpvB9acp3+wF4cOOu1aViFfF2Dtc9hr2I4TePwrUJujoWXPV1LLhfu3FWZtt1JZVb12Fgz+/hon6R+buverz6wFmWZvo5LxfKyA4Yy63b9tW2dTFHBrdr/LIqEiTsuG5eMu91mDaevx+v3ldhXcMXRaODxpVMfP72AACAU7LNzuxCjmCc8ODw9dX0TdmV/kmdrsuYtETbpLTiTZGMml3FzORfph0rlWV9j3flOLjZvb/f+j4ALmfOWBUj3l4u3i51DNMQw28BhWujVIU3q3Vh9v4ickr5276FpV1Ac6KUN5R0y3/33zpO1P3elZTLfDqLkWvZl9xd0HsJMOsJbyXR7XJvAdFp5MrJ7NWjhWCFDVZxHAvLkWDnCsaP1CX4uWle19OWdeTbWvax7n/3t043FKE0yOky3J/K/u2G11f596ztAQAAy5BMiLvH9yf45Te511feJZOj/6ZO12NEWqKWkFa8OJ/5Eu3MYEh3vkpm0H6Yie6n6lNJtQfAmTlWnYV4m+Yax5AYfusoXBtBL9DCnsiZrWZ7qPqqk+FgdoU7e7V03Ze0W3Wh2F5O9rZQ0i26TtTSX619VTRj2ca8hfXQAjYNNqdf/Ru2KyvM7iDbc5Btk8Ful11oaYqvUWCKqqJ2bZOO76uWZvmXxgU/eVkN8luhy8uraXXfy7Sd1W+dLN/Z+R+qmdvlvrll6mC/PmN7AADAUqIMxBB9UOY+VlKn6zYqLeGkpRUvLGS+NHnZTESFdKdmYiUtuZXtiWuNVC0N/INhL8p0jWnyJMsYc3iAxzdvrDoP8TbNFY4hMfzmUbh2SrkNtZr04tVqtvUFIhf0Ri5ofxHvZdzojAsnfXQhVOTv+Iuj3+tCnfaF0yfT9ZCo49etqg02UMgWtiszxdumbr9taaDaVYV6jVpdmflSl0a11lmU39yyZPq4mlfSstpkWtt3gP9zjkCVuD0AAOAimrUVxkudLskZacXL2Zu1z2hKprX3AaVmSFdrU0gCqJR1r2uSaEuDtWTWmn0QhXSqTPe1s9nWeiBtB+CiseoE4m2ay2w/MfweULh2wsmOCkXdLLPVVFFO1mqSVjVQGcmOph1O2i9ahW2hUOd0k9CYLfG3JfbtQrbjtudhu/JNT4eH9bLj9tVxE9TGOou+JpSpy4pp8BzXMeM02WbTfQwETUIBAFB1P7D14BPoQ789h7PSip3m39/7tR+netDZn6yRjJhNT1a1QGytEfl82BUh/RcyeErG9TVKyqLq+0j7PPL96zZqSrSbMQF3j9h4afcQb5dADL8PFK6dMKajwrh5YbNwyBcaNftV80EhyzfdhUry2f6l1TrtF1NoiX27kE3bnn9tlNyH7bIvGNAL6Hjw112j2WpvU8r+JpTJy7qIumBv7PYAAIDLSO1AOXW6FOelFZenNQ98OuvkW+Sy3LxFLSE8TbOGLkgkvRSnmPKdZtxkGvez9nmkuySTeRW2Zkn1PYBjl44HQ4i3aZZeBjH8flC4NqjuqDCVf0NI3I65CgpVU0NfFTP+3RfynFeo4wrZDr6UujRFT53O6gI6HrxmJ42Z2fjHBPGFGdW22wxc8V3L0cEb6hByKZ3Nd0duDwAAjy83O/sUPB5clw6Dv6WSTNOYW69kaKoUhO/AOXW6c5yfVjw24/7WB5suV5YVB3M4I01Td8x9rOovqbleh8PObHwOV1wjjQcsK/VavUas6kO8veljSAy/KxSuDao7+UtVv5Fyb6rKa775obvA8ry6QI9+rwrfzpZt6iqzXZ0QyoUUX0Rdw9FF7Nc5rKuWRbkqo+G3DinLuoRoe7auet+o7QEAAIsaqtEeOnDukDrddOenFRdTbutMmWacBtJYIdNVSua1+jSTOq14S31LAbficrHqNOJtmsW2nxh+dyhcG0sujKOCKS/Ucuo46UK11FIuLvlHzk57foZCG/nXfnC/y3/sBbFwP1++Rt3gdvXy66ybo1sz/AKG85Z1CbnZuGqyVQ3C6S+UAAAAc4kyUANph9A8KKSpUqebSWpacQmaKVsV1fpopuxUux5Zp2qtWv0ER5IyiT7dG6UdAVw5VjUQb9MsvP3E8LtE4doJ8Rs0+t6UEWo5dZ509RsptSDKn9RxoY3/bH8P/bHN1c9XXdqsEcDPsy7d9jXmpgnrrxfcvm5Cebz95y9rFhLYhsJJo4bhie0BAADLOpn+it8YF92sU6cbpSctcXKZYjitOLe9WftMWVaczpSp7NVUqbW+bkTqB4+6AaM2QdJUvuPsrPAvkAKgFo1VE6Wuy6LbcAfxdrntJ4bfrY+J/vZ//oOrDXOR8/NDN1177zvt8FFkbnyTfeS7g/teHWRemfvNfGRF/FvkUHxkdpz8I7fzyj8aSz76Pfvom9Wx3Yec6B9ZXnzsDgdZo8hBfovWvTnPaLsyWZ+O5R12hUzfWtegWq7d7sztg979mbqserrhQzUwXti3+ptbsO6no3WIj7MbRp0fAAA8I58O6EsnDNjl7l47NG18X5Y0TJR40DSDv7dLil/GjKVON+BkWqK5zKS04klj93e0LlO2URyKeD139bSSnmzs08bmSfpNx20krGTfyDFO2tfA3Uu4VifGuJCXHXNtEW8TXPMYRvOcst8EMfz6nrtwbWBoXIzRRd43ZIMFMfGFJ8PQRdT5+5C6kKt/0MDjRm+QaePlynhaUJY1vusPKs392LpQj6Qsq94vyYVr7X3rhq5gGwek09sDAADGaN5fB4b2TfxU+qvvAWDqdL1GpCXOTivOY/S+1uFofbq3sx460pMh494zkCkD+iXFqjjvd/w78bYaLhFvrZmPITH8vj1l4dqYk/a48OXwsSvyVmGQFg7JBdEo6e0WL/NUwU7X74O0NNquW2u7tPBKrqDh1evaLlkHmbYqxXajdYmDyagLb+qy6gAxHB9PjGdL66N9I8esexdHgY5AAgDALJIze1ZVE6GRefHpBjdGt9TpeoxKS5yXVpzFqYxSPPTt765tsC0k3CgNfWk72T/dEwBomBqrdHw/3nF+hXh7wXgbzHgMieF37UX/IztwtL/7X//Qfbq8v/h7/919ApZQmu1qZYpSQtCZrzoGAAAAAADPgRcaAF54s0xm8iVf1QoAAAAAAB4GhWuAs9/6t7LkhrI1AAAAAAAwBoVrgLU3e/fW4iz/YihbAwAAAAAAY1C4Bqj93lRla7nZ0NcaAAAAAAAYiRcaAAAAAAAAAIkmF64BAAAAAAAAqNAsFAAAAAAAAEhE4RoAAAAAAACQ6GGbhb68vLhPAAAAAAAAwHRjis3ocw0AAAAAAABIRLNQAAAAAAAAIBGFawAAAAAAAEAiCtcAAAAAAACARJP7XPu7//RH9+ny/uKf/cl9AgAAAAAAAK6PmmsAAAAAAABAIgrXAAAAAAAAgEQUrgEAAAAAAACJKFwDAAAAAAAAElG4dmf26xfz8iLDeu++Qb/SbFfV/mJ3AQAAAACAJTxl4Vq5XVUFVHZYme2Ugpdya1Zh2hez2pbuB2AB5d5s16vGOfeyWpm1nLTlg5965XZt1uut2XOJAcDdK/dbs5b7V7iX+ftZT5APDxPHDKut6ZpLuV+blXvIVg1yP33S+8pe0hKn9pfsscnHKZBEiS6jkV6x+1vSK24UAN3mjlXE28tbevuJ4ffhKQvX3hulEqUptn0n6LH9tmiMW5bv7hMwLw3SL6u1KSQYNs5PDX6FBPCv48/bu1Nuzddib/b7whYkAgDu1V4S+i+SQC8kk9G6a7nE/Py18aua6zZT0LyByr1V7is2o+G+egZyT92ezOHpcdKMVMJx0gfPLvOmU2ZZZrLqB5fhXMvcARybO1YRby/vAttPDL8bz90s1J845d58O3W+Knti64fc5Ln9BqdoAZEtMX/ggqAlyLn21QfArDC7w4f5+KiGw2FnirwKeQ8rezWv7mOe+08O5xQA3I/y3byXkt7Kd3L/qu9lH3IvC7cyievtNH++i8btHA6mcNNn+ReXEahoJqKwN4hM5nOopwnL1IzGhMzCXd93JOP3tXownA8lXhOPk3wpGTr34FmnlWkOh4P99+NQhHT2+uzcJfB4Zo9VxNuLm337jxDD78lzF669bszGnqNae+30CVN+c9Ui5cSmbA1LCueanGm7w6YOiiLLcrPR4C3fR18/GNlud1PYxRsPALgv2UYS6ZJQ3+X6TLMm97Ldm0u4i/1+YsJ9v60zNF+iGYcHoZpPODTvIY1l7s32Gbr28PspK1yat0ficSq3W5dplPF0WvvZ0XnKd5asBz2pAJElYhXx9rIusf3E8Lvy9C80CCXAcqINhxm5KKqoMlxqDMwgNF2W6NiqtwUAwA3Ym7Xtk+WMJ/LZl/rh0fu7e6g0Tsgg5BuziXID4eFUX0ZEMgv++3L/bdIyrytlf/vaBpkp3s54IBcfp4bSfPNNlfoePIfvS61YATyBcdfqxWMV8XaCWzmGxPB78/SFa/UJo/072Q/dfOHbqVLjSFeHgr4987CqXXTjxQlyYY27IH3b6Hra4Q4Vq3biOo4tjS51nevpmiXUsl5bnffpbQovjfBVSMvixIsgpq73PCYdI19N+Wi/OLLvqm3s+T0E6tNvL331jxzK0kyPY+P3ZThOQ1WvB7drqfPNj1fvqzHn1PnbAwC4eRLHQ20BScfFwsOp19fejMjYe+z0tIyMkpz2m1+okdDKEJ8jy7of+fV9L3s71KJ4f+emC3hzxqpFEW97LX0MieH3h8I1k2v5mrUfeLGBL7FvtzPvVhUMdHUoWMrfWnDWn/nXApgqKMS/2wIMLT0fvCp1uVqwsdfrV68S+b+urRaAVB0qDhYmfFvLeuk6u791urA8t16FzltGCPOut6kR9DTI6DjuT2X/dsPra/zLmeudJOEYnSi5r5tyyvRdnfj5AtronOuTffHVcmW/ayeSo7d/2r4My5EbmL9xttVPZXIT1wSfuqwjg+dbhxHn1HnbAwC4rPfq/qEGMidt4eVSRw88JfPi7iP9GQX5TZZV0T5q3Mcuk9MyqWm/BUiG+KttcVE19TmLZGBd443WNo+h+6n6xEvAAG/mWDUK8XZeCx9DYvhdonBN5BvXDrnvxQahxD43mxHFxnpBVydw1bGh7RDQDgez871Bamm8L5kPNFD46qeZKWTaejotNNgPFrKE5Urwsx3ga2eEMtjp7TVZmqL3DZPyW6FLlgvYd4Io04YgKhfTu25PUXV0WM/7wxzcNpXF11CYUnWmWLX7rr4ozJubRoc4Rpy33mnSjlFUENtRzTE8vRBdwaf0ETgU0g3INuYtrIcWsOkN5HQtvsn7Mqom3LVNOr6vLtzbgekS51uHUefUGdsDALiw8NBJb40DN4CGuqXBcRyPMo9D9GGM+zgkKS0zKV2xHJ8hzorN6TTHAK0ZEtZZ9ke8zbqtIdM1psmTpIPGHB7g8c0bq0Yh3s5s2WNIDL9PFK6pkCHvrnEUarqMKRSJq87Kxa4dG9YXlASAjQQAf9HLyd6o2VN+C4UnOu0mNI7W6TTg1J0RHgnLrdpkh0ktDTy7at0HavTYaeM3isRPD7Kqk8TdptXRocg2/qI/UfOoyyzrPdEZxyjcjKIbVKW++VhHv9eFOmNvaJmuh0RAv25VbbCBQrakfZmZL37Eo3UW4ZyU6Ts7MF3ofEuWuD0AgAvb128fk0zU0MOVWNz58tADz+lP589wRrpiEXv3RjjNnE5sS6SZVm1WpYNtYrXWDJ5sQ3EwH81cmRXSNHK//9qxYdryoj8dAOAysYp4u6TZt58YfrcoXLPqDHlZ+CDi+RcZZKYYEYlOdmwo6iZ/zcKok9NGnSK21QWAfW2y61pXfe2lbYFeUmzQarzu40RzrPdU5xwjXZlqklbVXluzT2T+91ZhWyjUqbdnDPsUxz6FaReyHfcnkLov4yaojXUW9b5qNqG87vk2rC7sHb89AIAh2jVE3adNNfia9kO/9duvo5r6oztqrl8uVd+Pr++sdEWnc/Z3abYuJ5TL/XDarU6bOcn6ldVQK82++Nrdj1FUE6IsXGbOZuyqddMaE2EqbfblPgKPYf7YuATi7ZBbO4bE8HtG4ZrTmyH3tV9GZsbHdGwoV3yorRMXPIyatkeYVvtmc6XV7cE/sDivvbRcnNutWUel4qvVV1c1d7q09ZZg3zGeHUY8njjnGEkEcoVGzVqOPtBnuZxHLkA1ptNAp/9qnwH2iyn0KUy7kE2Wv66b4arkc6C3KWV/E8rLnW8p6oK9sdsDALgcfRLu7xFZ8Tb+QYtPk0kEP/XAs3nvXtZ56Yp5lVuXJjtq/jNGZja+qwY3HA47U9gV13SH3N/9gYvkOxlPFuYftFYZO5mbpJ21i5Pp6wE8j6VjFfF2+XWbcxnE8PtG4VoQZcjDiw3qkuNxmfG6Y8Nr8qXV7cEb6nRxSPUWl5VZF4XE2+55nyOeXzx48Xp3jVcNpw7A+cfIv/UlbpteBfrMNjX01Wvj330hz3mFOq6QLTQPLk3RU0+3e99070ud78bfNMONVES17YaqgnctRwcv9Xw7R2fz3ZHbAwBoy80uSqxXg2v+P/hbB30o4xL3WXEwh9HxuE6T9dealszWmNlJRqi6S71qdzhnWiLtl7q/p7W2GEMzV5soc6X993Q9x6z6S2qul2bq6i5OrpMeAJaVeq1eKFYRb0e4pWNIDL93FK5F6jbH7sUGkzPjdaeAVyMXRnxRdA3jA2tEg7OrnnZ88R2Mb9qebNJ6dwU6N5wsWj//GIWqxuEFGHstwxEuaMp5ZNfi6Peq8O1scfPgro4lU84Bv85hXfX0d9WAw28dljrfzhVtz9bdQUZtDwBgOWXdsbJNS0y5P4Q0mYbx01F86GFbeMHQLG4g7eeFB0qlKVyTnsbgayyUhVm57zoqMXQKL/+SeXe+Eb1Xna64aL9MwJ1YLFYRby9mtu0nht89Ctdi+cYVElUnXXj1cEpmvKvQw4sCVnyS+hpRg9P2OGfaMULzOg3OUbXRcy293oMSjpEVqhqXVTt+HwjDeSL/2g/ud/mPDasjmxanOm9f+nX2x3r4BQxXPW6j5GbjSnyrGoTTXygBAJiRZvRWPl1VpSWmCGmygX52GhmvgftTaFYU7tszSU1XPLKQWazTGQAWjlXE2wvE2ytv/6UQwyehcK0hfrGB7zNqWrXMuvZb/1sSQw2a1kma+Su0b9p9/9s66ml9balldFcHHQhungSWrnL7S6137JxjVKnPEy2I8k8k4kIb/9n+Hvpjm6ufr/oJgkZ1P89z92VYfw2ie39T6g6k1zhuR3rOKa9Rw/DE9gAAlrQ3a5/Rk8za1IyeZhT9/frUvfTkPT6aV3zfHqXnvnNymWI4XTETyUR31SAPg9/vegzcd2MPRb3+EzKrcv/1byjMCt+3MQC1XKwi3l4k3orZt58YfvcoXGup3y7iTK1tFGq/SWhbawGdP42V60jQN69sn6Staev20FqLbmWrgsZza4hq3RVf16G0PqZtrNertDec+JpK9m2qjU3SAN69PEsu3mpKLbBxI0mwlP9XFl7vTuccIyecJ/u9BE0dtxW8/XaH36c0CZV9+qJvZ9F97dvpO35/2z9aBb/n7kvZgGpu2pTSnWvhu5ZrHDfv1DnlRTUMC/saatG3PQCAhZRm6+8FmiE4TH37mcwhyiid7KajdY+v7sEVvS+F2hyyLqOfnY5Oy1TLTElXXJvtV1eOk+6v5v20SoP69T/ab5J5XGuaoTGRzEO7E4mO+9s1uogAbllyrNKY4poJrnw/4R7x9qLxdpFjmIYYfiM+Jvrb//fz1Ya57HLzoZtu8p37Jnb4kIuk+l2GrDi475sG53EoPuT0C/PoGrLOZYtd3j9tln/s5Hf7uXP63UcerbtcCR9ZpkP8nczDjV2pt7dvlSoy7zAPGXS+4e9qOd3zaO5PPzT3a8p6jzGwbeccI6u1XVkh38RO/T6kta87h6zneJ23L8N57abtOf2dJc+3ofFa+9YNXdfqoajOy2o4tT0AgLk14/CJofPGEN0Th28ctVP3eE1PuVHHGXHfOTtdcQE+DdmRJhl1nLr2m59n3zAp/QM8maRYFecTmr8Tb+vhYvF25mM4iBh+86i5diR6c+KYEvsu2cYcPg5mV7T7Jsvk79zsDof+KrraLv6wM3ljwszkhfvefdOtmne9XC25rkqvM/li3Dz66Lz9q3yFzlf+yfLCLvOtd6ayP99kunh7ZB80a3Etud49zjlGVt00VB1Xmz71+xC3r+26taaSv22npLq/OlfvvH1Zd3YpTtbavMJxs8acU5VsEz2xOrk9AIC5hW4EUmlXBfbDhG46/D1e7sONsJ98bxpx3zk7XXFd2UbWT9Meus9a90pd/2JX/X60BfmmY5vraT4Sas4ATyMpVr3WrWVkvDjCEm+vEG9nPoapiOG34UVL2NznUf7uP/3Rfbq8v/hnf3KfAGAMrR6/MloTOpv0CnIAAAAAAMah5hqAxxXeFpRN6PMOAAAAAIDxKFwD8LDqV4nTJBQAAAAAsAwK1wA8KH3LUPVpWp93AAAAAACMR+EagMcUOmVNfDEJAAAAAAAjTH6hAQAAAAAAAIAKNdcAAAAAAACARBSuAQAAAAAAAIkoXAMAAAAAAAASPWyfay8vL+4TAAAAAAAAMN2YYjNeaAAAAAAAAAAkolkoAAAAAAAAkIjCNQAAAAAAACDR5Gahf/ev/pH7dHl/8W//m/sEAAAAAAAAXB811wAAAAAAAIBEFK4BAAAAAAAAiShcAwAAAAAAABJRuAYAAAAAAAAkonANAAAAAAAASETh2g3Zr1/My4sM6737Bv1Ks11V++u6u+tW1gMAAAAAAFzDUxauldtVVYhlh5XZTikUKbdmFaZ9Matt6X6oldu1Wa+3Zn/8E27EXRyjcm+261XjfHtZrcxaTtjywc8triEAeGx7ub/V97atGRvup0xX7tdm5R6AVYPcU5/s3lJKWmItaYd6H8ggf6/WkpZw43jhIe/A0PsgURImemwaaRa7v4+XA6Bp6VhFvF3eUttPDL8vT1m49t4omShNsZ0QZLZFY9yyfHefnHJrvhZ7s98XthAEN+gOjpEG6JfV2hQSkRvnpga+QoL31/Hn7N3hGgKAxyZxfpuS4xg9XVWr3GYKmjdRub/KvUUfVD387UX2gWaUJC2xdzshyzIZ9CfdD1vzba6EhByXld2nVZrFLqf6wWU4ZR3s3wCaLhCriLcLW2r7ieH36LmbhfoTp9yPOzltkNEPuclz+82x7NW8uo957j89KC0A0lLtCU9AbsKtHyM5z776KJwVZnf4MB8f1XA47EyRV+HuYQ0dn3s95wAAjmQYvlYPKvPexFSX8dPp0/fC3iQyk+8O4R76IffQ6haqT+gnZBbu8N5j94HNGGemcPvgcDjI4NITu4350pOcyHd1uqM97I52vdaqcA+e8505yDh2OTr+oQjp7PVj566BJLPHqiPE26UtdQyJ4ffpuQvXXjdmY08wrb12+oQpv7lqkRJk+sNMbnbh5O0543Flt32Mwnmm63nYuMBcybLcbDTAyvePe3ZxDQHAw9pvq4xIVrg02EhjpwsPQuVuIvfLxn1E7qG7N5dZkAzFtqNrj4cgmdMqHySZMskkbTrupZmkZee4w5bbrcs0yr6VXFtjntlGMoDuYMnxe9TdDSS5RKwi3i5rqe0nht+tp3+hQSiN3+/didVHLoqqWHpiyT8wTWi2nGWhBhcAALdlb9a2T5YpT+T902/JMLxNeUg0frrwgKovUyiZBf99uf/mHmbdg7H7u5RMXDVGVryZzRy5r16l+eabjfU9eA7fl6bdkwrwmMZdq8vHKuJtumseQ2L4PXv6wrX6hNE+nuyHbr7w7VQJvl4QrjPDqsS5TauGNjsLHN9RoG/3XE873FmiX5dVVdJcbs06TOu+C2S9tjrvZoeJvv11LLwQwm9gWZx4ycPU9U5R7dd4/rbjf/dr0/Ax0jbs9X7qHrqPrZqyHt1ebWN6UZZmegwbv6/DcRyqdi3nTHVs2+eLWup8PD4+Y86587cHALCk8IQ8lwzHhAzDlOnCA6rX195M4dj77PT0joxi0xBxOqA7LbWcd9fvT2byvjZDC8iyvseBchzcary/c+MFvDljVRfi7fKWOYbE8HtG4ZrJtXzN2ksw6Ttl9q7kLcu/9F48p2kpeHXBx8uxBRRaMj54xWmBgxZc7KsLTs5y7WxQv/edJQ4WFnzTDvKLqNBDpgvLc+tV6LxlhDBvjQOuMDCeuQYQHcf9qaoOFqvh9TX+5cz1HuNdC0zagVTWWzv+9wF6JLuta7+f6v0QuPVvfVuZaT2yL75KrhwX24ml/WOEafs6LEduXr5Kc1v9RCZvtetf8nzsMOKcO297AACLKvVlNRqBq6Yno02aTjIv7l7Sn1GQ3+SeUnmXzJH72GVyekcydjYN0ZypT0tdpA+h0BLjVTJ19sOV6X6qPh29BAx4WjPHqjbi7QXi7ULHkBh+1yhcE/mm7qyv88UGEmiqzHpuNsl1MzUI+KqlWeiY8OPjYHaFFgrIhTRwwWmQsLEucx3c2w4NdR7avlvHKE3R+wZJ+a3QJUug9J3j2/bb1a9aqvEua5AXVSeG9bw/zKGotrcsvobCkkw7O9RxfNCVdXpz0+gQx+Lz1nucfaGdMMadSEZt3ve+zfoI9jjblXXHp72u4nVjv+uqvTjbemQb8+b2u56TWrtrTC2/yfs6+xL6c/OFx011VeF2ofKi52OHUefcGdsDAFiWf9t6Vmwk+o83bTr/xP8EfSDjPg5JSu+4dIBNT/m0QLinF/ah1JLKOrdnXiUNYd82F9XqqGqYuwdNPXRbtPaHH9Y2DdI1RZTpGtPkSdZtzOEBHt+8saqNeCsWj7fLHENi+H2jcE2FTHlp9h2la6G2S19b5DHKb6FwRINA3TGhBIWNBhPf4WGHULhXtX0Pk1oaVHbVeg3U2LHThreW6J/RhZ5tJChJMNq0OjEU2cYH1xM1i7rMst5jaCGNBt2wcTLvN+Nj69gqrPVx1vbx8crqMXLHJzxN6DLPeqhsU91c/Jyq2mADhWxJ+zozX/yIXdsVzlmZPq7mtfT5mCxxewAAy/IPmDSzNOUhZep0olnLYWHhvlil8TQdUC9d7jl6T/eJgUt1Cq0Zy9Xavm1OH/zVNUK0hrlrMWH/7iLjSEbMD3tXI/0lPCSuhX6IZXlfOzZMl9WfFgAwe6wi3l483i6y/cTwu0ThmlVnysvCtTMP/IsMMlMMd7Y2aEqHh22NQp/Oa7du2tpXgGML9JKu+7qd9VRzrPcYtjDnaP6Xr8I693rYJzj2CUy7kK3d9DR9X8dNUNuVvepzttmE8rrn47C6MHj89gAATvGdO8eDT6QP/abK0DlzLjF6fPhNne7y6vtLf7+89f1WMjsnkwPn7G9PMmI2HVHV5rC1P+TzYVe/va7qtLyW7+R3X6PcD/rQMNQEkWnaTa1kGb5GSVlUfR9VtSWqdWvUsJAE0S0fR2C6Oa7VORFv1eXj7RKI4feIwjWnN1Pua8CcmSEf0+FhnzCtllLbk/148NfWeQVJpdF+59aNqqRfXbXb6dLWe1/15dU1THz0EDqQHCm0iT96yqE3nKqatJYa9cTxXlPXo0mfwLQL2bQvgbqZrko+R3qbUvY3obzc+ZiiLtgbuz0AgOWUW5eOiBLwY6RO553z0G6qUWm86H57kXWTdOtbVAPey/Ko6wmfxo0cJVnkC60J8lHnvo5qMdgMnS7LTVvVltBJc9vNRsrxA57FnPGAeOtcON4usgxi+F2icC2IMuXhxQZ1Kf6tZMjjKp7x4A11qDikekPLyqyLQq7T7nmfI55fPHjxeneNVw0LF9RIsKpiVWkKW1pfFRTpv7aLMDkDzqm9eB5XyBaaD8s6tiOj073v+s6RzGz8NsUBOjShHO5nsGs5Onip5+M5QvXmhO0BAHTJzS5+Em4H1wXA4G+ptf9Tp5PM1pgQLxmh6k41R4fRdafW80nd3+PUHWxPED1c7MpIVv0lNdfrcNg1utm4RpoAWFbqtbpErCLeprmlYzgOMfx2UbgWqdscuxcb3FqGXE76+ITvGg4p66k1kGxQ7bqwDqHPsGST1rsriLlh6WLz/dY9tZHgY6OlLyjSarnacf8yTRkniZsPd3UqmXKOhGC71/Ioq276WQfiI0udj+eKtmfrqveN2h4AwLzCQw7/0Ko1hGrOReiw2X6VOl1k6IFc6DB6FnX3D7cgZLok/bLwI8mR6rTFRftlAu7EbLGKeHs1c24/Mfy+UbgWi2ou6YsN/BtT5siQh6aBXQUiJ5wz7Rih+ZwWrEVVQs+19HrPrdoP+tRmZwvS6kKig+yXdsf9t+W8fS3ntzvBq31QN6EMBc6R2z+uudm4EuHqzTfD2wMAeBRRxmvgHhWaFc39wGXovhge2Mp9dMkMisy7mvu7VpjoNG9m94SQca/TGgCuHKtm8ezxdqHtJ4bfNQrXGuIXG/h+o6ZWke1Wl0L3vEFx3/8mjnpaV6NuId1VPQcCl9dTsn6p9Z6Hr2asBatRx403p356oBHd3y7O3deNppT7usZmVxC9iePac855oSNTXccT2wMAWMipGs6+RnpWuE6bP6q+XVKnE3UrhJ70Vvymuak3hZ57z8llilCDeul7kaTlqrt0X/cR9QMnXZGxq1Ju/Qu/svFv3JZ7sO9wOyt838YA1OyxingbXCreLrL9xPC7RuFaS/12EefMFxkEoVaclqOtos7otYbcyla5DV+1xX2BfV2HkvhYud+adcfrdcfwNZHsm1LjeesFpfPsW7FQsq4FMm4kCYTy/8rC6z2v+ulDWaxDtWc7+I7617p/+nbGXGSfyzL9shpL88fD/tEq9D13X4fgrE0pTzShvOZxPXXOeaEjU1nH9Xw1UAEAN66d3opuUnpvWq3cPWHgTXNHRqd3qmU23+itLyGSNITvfmPxDEpde1sf3K78PV1JOkJfGlWtSjMdsV/rG+L0QW877dFcf93WRo8Pknlc6zIaN2KZh33pkUsHyL5+u3q/GsCNSY5Vek36PErrzY+X9uzxdpFjSAy/ax8T/e3/839cbZjLLjcfuukm37lvYocPOZ+r32XIioP7vql/HvX0Rz/t8g85LcO8G0OWf+zkd/u5c712H3m0XnKWf2SZDvF3Mg83dmVgXRpk3mEeMuh8w9/Vcrrn0dxXfmjus5T1HuP0th2Kar3HH6N4e/x6unlEQ3P7zlmPLq1j0TlkPcs6b1+Hc9pN23PqO0uej0PjxceoHrqu07Df7XBqewAAF+fTPVkh0X2CU9Mdiv70lg6a5nKjjjPi3nNqmTJko9IBc+he33o4Tkc00wDdg67/0f72x6JvmHpsgWeSFKvivMKEWEa8XcYix5AYfq+ouXYkenuilhzPWUqrfZoddq6zfC8zeeG+d990y20/YLvC94mmJfhVKX4mX4ybRx+d984UVVUfW8KtZde+E/+33pnKvnqT6eLtyfJWVdMl13tme/9CA60efZD1qga5TuzfO4lyytbws5+W4I6F3V/xfhTyd/XCib5XIp+3r/ONfxOpOFlj81rHdcw5V8k20dOquWqgAgBuX7YxB71vS7qmEfqT708j7j1+meG+6Mk6yLh6z9R+bS9D1rdxj/Y0HVGl7dqrku90/QqbRm1MEk1j++V13wb5pmM5MpVscyHz/DhsjqcBUEmKVa/aGKMi47mGN9fz7PF2kWNIDL9XL1rC5j6P8nf/6h+5T5f3F//2v7lPwNy0OaZWfZWA3PtKez9OZgoJUNSOvXVlqDqdFXJD4YABAAAAABZAzTVAhTehDCjfBzvRx40JbwrKxnfcCQAAAADARBSuASrqPDN06B/TDiS/+k4paWJ4D/ZbjhcAAAAAYHk0CwWcozfO+Dbr2peY/ULYdvxT+w7A5fkmvHLIaBIKAAAAAFjQ5MI14JHpa5O15tp7XKAmtFPHfLMxX9qdVeI27dfmZa1Fa0N96AEAAAAAcD4K1wAAAAAAAIBE9LkGAAAAAAAAJKJwDQAAAAAAAEhE4RoAAAAAAACQ6GH7XHt5eXGfAAAAAAAAgOnGFJvxQgMAAAAAAAAgEc1CAQAAAAAAgEQUrgEAAAAAAACJJjcL/fcvpft0ef/yI3OfAAAAAAAAgOuj5hoAAAAAAACQiMI1AAAAAAAAIBGFawAAAAAAAEAiCtcAAAAAAACARBSuAQAAAAAAAIkoXLu60mxXL+bl5cWs9+4rAAAAAAAA3IWnLFzbr6vCrJdTpVn7dTXey8psS/edU27XZr3emn3r+5tS7s12vTIruw1uWK3Mers35S2v9wzu4vgAAPCASkl/rCW9EdIeLv2xknTX0G05ZbpS0mor95CyGnT8Z7v/l5Jk9fttbYZTt6Xss23nfl6P2Gmp+5vjBPRb6vogFl/alFg8hDh9tz4m+nfmcLVhLrvcfOimm3znvumxy6vxTPZRxIs/FB+Z/X7EPE46fBRZNa+zZxU5hHXvGbJClvygZj0+AABgHEnT5FkjvZFlmQz+71Z6KkiZrk4/dQ/ZrOmqW3XYFR95Yz/kH/2bvWuN2zH07rTU/c1xAvotdX2kxFSVMh3XuJoWi4cQp+8ZzUJTZK/m1X3Mc//phpRb89XXyssKszt8aCGqHQ6HnZGgWf32qIaOz97VRlxtB5/YAACAafSJfWEfcWem2B1cuuMgg0uD7DbmS0cSJGU6O429kWeSz6imsYOkc6pkjtYgmFBz4O7SB1UNidW6sLUKJPN7Wvlu3svMZPku7NvmPhOyH3wSMpa6v2c/TsADWer6sPMlFl9IQiweQpy+b7LjJumqUXapYS5n11ybVV1SPFeJ8KHwTxxSS8wfmD+mj1xzDwCAS0tNM6VMF9VQ70w7Rb9nY2d6Z+mDkJbNJK0nObB6Pyam/YZq/afu7yWOE/Aolro+iMUXNXssHhLtT+L0baLm2gN69x2qZVmowQUAADDN3qxtvyunnliXZrutxsiKN7ORlPg4adOV31y/P1lhNrn9qinbhO/L/bdq3Lswdn8b85rlVQ0DW7Ng9A7vl32pa0W0pO7vxz1OwJBx1/Ey1wexeB5XjMVDiNM3j8K1JMNv+Kw6IKx+7xu6pqtUVUvrcd0LCNyvY7z6C7sszXv1aQLtQHFcZ4bl1q3nULXdcuteqHD8Uogpy6r4/e7mJfOu93M8/+PjE9Y1fFE0XvSwkonP3x4AAJ7RuyY5RGbyrrZGvdKmCw8RX19lym5j00Jj0gdtXR1Nr0Z2ND2XbLMzu75c1pmyrPloNnV/z3mcgEezzPVBLH6kWDyEOH2bKFybWdzmWgNU5k9GT/7W7zpP3nctuGkHhNLsi7XMs7c07kj2JXfz1ze9rHsKqrpooZSu/74KsG5d9ftyX9jgFce1sBwJfu5hx5G6NDxvtdGftqwj39bmZeX3s5Lphq52DRit/W7/dsPrq/x71vYAAPCk9nv3dP9VEuH2wzhJ00ni3t3v25mLWCb3/Yr2X+M+dhmRPqhVD++qdF5zpqX8bR+O3mufrpJJrfrc0V3S3Oa0/T3zcQIeykLXB7H4/mPxEOL0zaNwbU7l1mxtaY/vBFI7gKz+3flqlq8b+11Xtct9UcgpHnckqNO5C6en48JO2ca8FW46+yplDT6nX6MbOjPM3EsQbAeW8fqXpvgaBaqoaupegvKx0nxzC83yL41gOXlZDfJbocvL65c19OxTr+oUUpbhD4Qs980tUwf79RnbAwDAsyrr1Ll5lXTHVh80RjUJqlrpx7Xw06bzNSxOkIzHmPv0qPSBE9IuLq128J0+a9olpLsKu873RGt/hHWW/RFvc/r+nvc4AY9lmeuDWHzfsXgIcfo+PHfhmn8bSd8w8YIMtZpybZscn4Jy4W+K6qQMTwa6aGGRBg8/rQaMN+NjxPuEouJsUwUnP6eqNthAIZstGNQPmSneNq323LoeO1k70ajVlZkvdWnU8XaV39yyZPq4mlfSstpk2vitKXIzaMwmSeL2AAAAe99erdb2bXP6sLCuhaC10teSWevpvyZxuuaT+4WFtIvma6q0Wr10SRdoussn2CQTdMtdR2jGVJtO6aDpXa39YR/uFgejvXP3Sd3fFz1OwJ1Z5PogFld/3ngsHkKcvk/UXLshtlDp6LzWarDVp3Kw3eMx+wTAluC3C9mO26I3CgY7r63c5O46jgv54iao7cpefU0oU5cV02A6vpPO8bLNpirYm7A9AADcP9+Bczz4jNTQbzHJ5Ni0R1WDwNY4kM+HnXvAqPPpfHCZOt3l1GmAnk6fRZ0mkozoySTbHPs7hTYDkvUrq6Gm3ZB8vWhfRcD9u9Z1fAqx+PZj8RDi9L167sI1CR5VFdKeYaBUuEtoj3xUSl6a7VZLm0Weu8Kb8UIngkm0BL9dyKZt0b821jF0ZqhPJVwpeXvwsbRRyNfblLK/CWXysi6iLtgbuz0AAEBkuXmLas17WR51V9FVMzxxuik1+s81ptPnOE10yXWbJjMb36WGGw6HnSnsild9FfU1pUrdptvdF8D1LXJ9EIvvIBYPIU7fK2quzUkCTxV3SlPYt1VWBUX6r+0iTC6UYqhjsEW5QraDf/Ig69jT5tKXkrcHr9nhoVz8fpviYBuaUObye2/o61yODt5Q54pLyevStcnbAwDAfcrNLkrIV4PrpmHwt3HqDpGnOZ5OO712H4dIRqFKTUzs2LtT3enzfJbd31NkkqHe2OZV1d/at0/9ADZ1f1/jOAGXlnodX+/6IBa33U4sHkKcvg8Urs1JTnLbuaLWTrNnqpyItqBIq9hqx/3LNGWcJNP+4NxniU7hmvRO1eaT4dDeiFAbb6/lUVbd9HOgpl7Ksi4h2p6ti1qjtgcAgCcUMl2S5pmS70mdzhuq3R466J5F3UXHIwv9A0uKZ//tKIWYvL8vd5yA+zPn9UEsfnzE6dtG4dqMqmaEWjttZwvS6kKigzns2h3335bQ9FQunuPL9JRcy5ysah/UTShDLbDIecu6hNxsXNXncv9N1nF4ewAAeGqvmUvs97+mvzNxnjRdlLkaSEeEpkNzPxQbSruEWu6yaQ/zeD91f1/5OAE3baHrg1hcechYPIQ4fSsoXJuNr6aqpciultNNqmuX6dXkw039xGJvOgrBT2o0pdzXTSj917FzlzULCRIdt5YgdISp63hiewAAeGrZq6nu7H1dTtQPqRqJ88TpQpqj763i8dvkpt64e9IHJ5cpQi33O00v1OvfzJCm7u9FjxNw5xa5PojF1r3H4iHE6dtG4dps6pLfstBXFUdvFfEd9a+3Zu9LfhdTvdXEL6uxtFJ+W/k3nLT6f4v7i/sq43SsprbtrqdvkYutmps2pXQXffiu5dxlnSM8mdFCQLdg3U/tdQgdYco62lcfi77tAQDgqdU1vu3Linw6QEnaYytpINttRjvtkTpdSEfoZCtJ+IepbPphtXL37YG3yR05lT5oLbP5tjbXwXS1srJY/+bx21JuV5ImXdv91Uz3yPrLb379j/Zb6v5e4jgBjyL5+tB44/OY22qcgFis3916LK7Ws/sYEqfv3MdE/84crjbMZZebD910k+/cNz12eTWeyT6KxuIPH3IS2t+as6i/12myrBqqv+sha8ysb161Q+HmcWp9rd2HnPeN5R0PWc+yZNqw/tV41TbE3+UyVrewX920zX3WlrKs0/uqMjRefIzqoXlMKmG/2+HU9gAA8My676/10Jf2SJzuUHxIfqBjfDdk/emVbiPSB6eWKUM2nECZTTONMjBE6zNqmr79lrq/Zz9OwANJuj7ivF7X78RiHW45Fg8dQ+L0faPm2pz2/oUG2lH/wRwO1SD72f69k0ihymI7f42sIDc7fVVvkZvM923myd8SaOw6+TeNNOm0up46rf6tJeZVqbnOKy902v43pNQdLIosN19ai286b1npMrN5k/1TLbQi65p3rGy2iZ52nNweAACemdxfG/d1T9Me1UudutMeidNlG3PQtFXuazk4yWmIEekDv8yudZVxdV0P3Rt5E7KNpktlG3WfNdZffpP1L3bV792HKXF/z36cgAeSdH28amOaioznGk5FiMW3HouHjiFx+r69aAmb+zzKv3+pqwpe2r/UMtWbpc0xtRmjXNC9r+j142SmkIv+6m8OxQllqAadFRLIOGAAAAAAAKCFmmtz0Y783cde5ftgJ/q4MeFNM1lnzTYAAAAAAAAK1+YSdb4YOvSPaWeQX32HgDQxvAf7LccLAAAAAAAMo1nojI7eWJK5AjftS8x+IWw7cNot3z7fhFcOGU1CAQAAAABAj8mFaximr6zVmmvvcYGa0A4I883GfGl3FIjbtF+bl7UWrQ31oQcAAAAAAJ4dhWsAAAAAAABAIvpcAwAAAAAAABJRuAYAAAAAAAAkonANAAAAAAAASPSwfa69vLy4TwAAAAAAAMB0Y4rNeKEBAAAAAAAAkIhmoQAAAAAAAEAiCtcAAAAAAACARBSuAQAAAAAAAIkoXAMAAAAAAAASUbgGAAAAAAAAJKJwDQAAAAAAAEhE4RoAAAAAAACQiMI1AAAAAAAAIBGFawAAAAAAAEAiCtcAAAAAAACARBSuAQAAAAAAAIkoXAMAAAAAAAASUbgGAAAAAAAAJDHm/wensUSdq+MovgAAAABJRU5ErkJggg==)
# this code takes the above function (remapDNBR) where the dNBR threshold has been applied to the image'
# and applies a color coded map to each threshold as shown in the image above'
def showDNBR(dnbr):
cmap = matplotlib.colors.ListedColormap(
[
"blue",
"teal",
"green",
"yellow",
"orange",
"red",
"purple",
]
)
plt.imshow(remapDNBR(dnbr), cmap=cmap)
showDNBR(dnbr_remapped)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[35], line 20
6 cmap = matplotlib.colors.ListedColormap(
7 [
8 "blue",
(...)
15 ]
16 )
17 plt.imshow(remapDNBR(dnbr), cmap=cmap)
---> 20 showDNBR(dnbr_remapped)
NameError: name 'dnbr_remapped' is not defined
Further Reading#
A Project for Monitoring Trends in Burn Severity | Fire Ecology | Full Text
The unequal vulnerability of communities of color to wildfire | PLOS ONE
Fire Frequency, Area Burned, and Severity: A Quantitative Approach to Defining a Normal Fire Year
Science, technology, and human factors in fire danger rating: the Canadian experience
Socio-economic effects of wildfires
Zhao, J. et al. (2020). Quantifying the Effects of Vegetation Restorations on the Soil Erosion Export and Nutrient Loss on the Loess Plateau. Front. Plant Sci. 11 https://doi.org/10.3389/fpls.2020.573126
Amanda K. Hohner, Charles C. Rhoades, Paul Wilkerson, Fernando L. Rosario-Ortiz (2019). Wildfires alter forest watersheds and threaten drinking water quality. Accounts of Chemical Research. 52: 1234-1244. https://doi.org/10.1021/acs.accounts.8b00670
Alan Buis (2021). The Climate Connections of a Record Fire Year in the U.S. West. Link
Ian P. Davies ,Ryan D. Haugo,James C. Robertson,Phillip S. Levin (2018). The unequal vulnerability of communities of color to wildfire.PLoS ONE 13(11): e0205825. https://doi.org/10.1371/journal.pone.0205825
Wildfires and forest ecosystems