diff --git a/.github/workflows/dev_ci_cd_conda.yml b/.github/workflows/dev_ci_cd_conda.yml index 2e59d6f..c8b81f6 100644 --- a/.github/workflows/dev_ci_cd_conda.yml +++ b/.github/workflows/dev_ci_cd_conda.yml @@ -33,4 +33,4 @@ jobs: GOOGLE_APPLICATION_USER: ${{ secrets.GOOGLE_APPLICATION_USER }} GOOGLE_APPLICATION_CREDENTIALS: ${{ secrets.GOOGLE_APPLICATION_CREDENTIALS }} run: | - pytest tests \ No newline at end of file + pytest tests diff --git a/README.md b/README.md index beea254..56ae15e 100644 --- a/README.md +++ b/README.md @@ -43,6 +43,7 @@ To run the module, 1. You need access to Google Earth Engine 2. Install + 3. If you want to use the ERA5 layer, you need to install the [Climate Data Store (CDS) Application Program Interface (API)](https://cds.climate.copernicus.eu/api-how-to) ### Interactive development diff --git a/city_metrix/layers/__init__.py b/city_metrix/layers/__init__.py index 7e5e19e..e701e52 100644 --- a/city_metrix/layers/__init__.py +++ b/city_metrix/layers/__init__.py @@ -19,4 +19,5 @@ from .alos_dsm import AlosDSM from .overture_buildings import OvertureBuildings from .nasa_dem import NasaDEM +from .era_5_hottest_day import Era5HottestDay from .impervious_surface import ImperviousSurface diff --git a/city_metrix/layers/era_5_hottest_day.py b/city_metrix/layers/era_5_hottest_day.py new file mode 100644 index 0000000..92c2982 --- /dev/null +++ b/city_metrix/layers/era_5_hottest_day.py @@ -0,0 +1,115 @@ +import ee +from timezonefinder import TimezoneFinder +from pytz import timezone +from datetime import datetime +import pytz +import cdsapi +import os +import xarray as xr + +from .layer import Layer + +class Era5HottestDay(Layer): + def __init__(self, start_date="2023-01-01", end_date="2024-01-01", **kwargs): + super().__init__(**kwargs) + self.start_date = start_date + self.end_date = end_date + + def get_data(self, bbox): + dataset = ee.ImageCollection("ECMWF/ERA5_LAND/HOURLY") + + # Function to find the city mean temperature of each hour + def hourly_mean_temperature(image): + hourly_mean = image.select('temperature_2m').reduceRegion( + reducer=ee.Reducer.mean(), + geometry=ee.Geometry.BBox(*bbox), + scale=11132, + bestEffort=True + ).values().get(0) + + return image.set('hourly_mean_temperature', hourly_mean) + + era5 = ee.ImageCollection(dataset + .filterBounds(ee.Geometry.BBox(*bbox)) + .filterDate(self.start_date, self.end_date) + .select('temperature_2m') + ) + + era5_hourly_mean = era5.map(hourly_mean_temperature) + + # Sort the collection based on the highest temperature and get the first image + highest_temperature_day = era5_hourly_mean.sort('hourly_mean_temperature', False).first() + highest_temperature_day = highest_temperature_day.get('system:index').getInfo() + + # system:index in format 20230101T00 + year = highest_temperature_day[0:4] + month = highest_temperature_day[4:6] + day = highest_temperature_day[6:8] + time = highest_temperature_day[-2:] + + min_lon, min_lat, max_lon, max_lat = bbox + center_lon = (min_lon + max_lon) / 2 + center_lat = (min_lat + max_lat) / 2 + + # Initialize TimezoneFinder + tf = TimezoneFinder() + # Find the timezone of the center point + tz_name = tf.timezone_at(lng=center_lon, lat=center_lat) + # Get the timezone object + local_tz = timezone(tz_name) + # Define the UTC time + utc_time = datetime.strptime(f'{year}-{month}-{day} {time}:00:00', "%Y-%m-%d %H:%M:%S") + + # Convert UTC time to local time + local_time = utc_time.replace(tzinfo=pytz.utc).astimezone(local_tz) + local_date = local_time.date() + + utc_times = [] + for i in range(0, 24): + local_time_hourly = local_tz.localize(datetime(local_date.year, local_date.month, local_date.day, i, 0)) + utc_time_hourly = local_time_hourly.astimezone(pytz.utc) + utc_times.append(utc_time_hourly) + + utc_dates = list(set([dt.date() for dt in utc_times])) + + dataarray_list = [] + c = cdsapi.Client() + for i in range(len(utc_dates)): + c.retrieve( + 'reanalysis-era5-single-levels', + { + 'product_type': 'reanalysis', + 'variable': [ + '10m_u_component_of_wind', '10m_v_component_of_wind', '2m_dewpoint_temperature', + '2m_temperature', 'clear_sky_direct_solar_radiation_at_surface', 'mean_surface_direct_short_wave_radiation_flux_clear_sky', + 'mean_surface_downward_long_wave_radiation_flux_clear_sky', 'sea_surface_temperature', 'total_precipitation', + ], + 'year': utc_dates[i].year, + 'month': utc_dates[i].month, + 'day': utc_dates[i].day, + 'time': ['00:00', '01:00', '02:00', '03:00', '04:00', '05:00', '06:00', '07:00', '08:00', '09:00', + '10:00', '11:00', '12:00', '13:00', '14:00', '15:00', '16:00', '17:00', '18:00', '19:00', + '20:00', '21:00', '22:00', '23:00'], + 'area': [max_lat, min_lon, min_lat, max_lon], + 'data_format': 'netcdf', + 'download_format': 'unarchived' + }, + f'download_{i}.nc') + + dataarray = xr.open_dataset(f'download_{i}.nc') + + # Subset times for the day + times = [valid_time.astype('datetime64[s]').astype(datetime).replace(tzinfo=pytz.UTC) for valid_time in dataarray['valid_time'].values] + indices = [i for i, value in enumerate(times) if value in utc_times] + subset_dataarray = dataarray.isel(valid_time=indices) + + dataarray_list.append(subset_dataarray) + + # Remove local file + os.remove(f'download_{i}.nc') + + data = xr.concat(dataarray_list, dim='valid_time') + # xarray.Dataset to xarray.DataArray + data = data.to_array() + + return data diff --git a/city_metrix/metrics/__init__.py b/city_metrix/metrics/__init__.py index 10b559c..d95cfa5 100644 --- a/city_metrix/metrics/__init__.py +++ b/city_metrix/metrics/__init__.py @@ -3,4 +3,5 @@ from .built_land_with_high_land_surface_temperature import built_land_with_high_land_surface_temperature from .mean_tree_cover import mean_tree_cover from .urban_open_space import urban_open_space -from .natural_areas import natural_areas \ No newline at end of file +from .natural_areas import natural_areas +from .era_5_met_preprocessing import era_5_met_preprocessing diff --git a/city_metrix/metrics/era_5_met_preprocessing.py b/city_metrix/metrics/era_5_met_preprocessing.py new file mode 100644 index 0000000..513aca5 --- /dev/null +++ b/city_metrix/metrics/era_5_met_preprocessing.py @@ -0,0 +1,68 @@ +from datetime import datetime +import pandas as pd +import numpy as np +from geopandas import GeoDataFrame, GeoSeries + +from city_metrix.layers import Era5HottestDay + + +def era_5_met_preprocessing(zones: GeoDataFrame) -> GeoSeries: + """ + Get ERA 5 data for the hottest day + :param zones: GeoDataFrame with geometries to collect zonal stats on + :return: Pandas Dataframe of data + """ + era_5_data = Era5HottestDay().get_data(zones.total_bounds) + + t2m_var = era_5_data.sel(variable='t2m').values + u10_var = era_5_data.sel(variable='u10').values + v10_var = era_5_data.sel(variable='v10').values + sst_var = era_5_data.sel(variable='sst').values + cdir_var = era_5_data.sel(variable='cdir').values + sw_var = era_5_data.sel(variable='msdrswrfcs').values + lw_var = era_5_data.sel(variable='msdwlwrfcs').values + d2m_var = era_5_data.sel(variable='d2m').values + time_var = era_5_data['valid_time'].values + lat_var = era_5_data['latitude'].values + lon_var = era_5_data['longitude'].values + + # temps go from K to C; global rad (cdir) goes from /hour to /second; wind speed from vectors (pythagorean) + # rh calculated from temp and dew point; vpd calculated from tepm and rh + times = [time.astype('datetime64[s]').astype(datetime) for time in time_var] + t2m_vals = (t2m_var[:]-273.15) + d2m_vals = (d2m_var[:]-273.15) + rh_vals = (100*(np.exp((17.625*d2m_vals)/(243.04+d2m_vals))/np.exp((17.625*t2m_vals)/(243.04+t2m_vals)))) + grad_vals = (cdir_var[:]/3600) + dir_vals = (sw_var[:]) + dif_vals = (lw_var[:]) + wtemp_vals = (sst_var[:]-273.15) + wind_vals = (np.sqrt(((np.square(u10_var[:]))+(np.square(v10_var[:]))))) + # calc vapor pressure deficit in hPa for future utci conversion. first, get svp in pascals and then get vpd + svp_vals = (0.61078*np.exp(t2m_vals/(t2m_vals+237.3)*17.2694)) + vpd_vals = ((svp_vals*(1-(rh_vals/100))))*10 + + # make lat/lon grid + latitudes = lat_var[:] + longitudes = lon_var[:] + latitudes_2d, longitudes_2d = np.meshgrid(latitudes, longitudes, indexing='ij') + latitudes_flat = latitudes_2d.flatten() + longitudes_flat = longitudes_2d.flatten() + + # create pandas dataframe + df = pd.DataFrame({ + 'time': np.repeat(times, len(latitudes_flat)), + 'lat': np.tile(latitudes_flat, len(times)), + 'lon': np.tile(longitudes_flat, len(times)), + 'temp': t2m_vals.flatten(), + 'rh': rh_vals.flatten(), + 'global_rad': grad_vals.flatten(), + 'direct_rad': dir_vals.flatten(), + 'diffuse_rad': dif_vals.flatten(), + 'water_temp': wtemp_vals.flatten(), + 'wind': wind_vals.flatten(), + 'vpd': vpd_vals.flatten() + }) + # round all numbers to two decimal places, which is the precision needed by the model + df = df.round(2) + + return df diff --git a/environment.yml b/environment.yml index 0d00195..0350113 100644 --- a/environment.yml +++ b/environment.yml @@ -21,7 +21,8 @@ dependencies: - geemap=0.32.0 - pip=23.3.1 - boto3=1.34.124 - - scikit-learn=1.5.1 + - cdsapi=0.7.3 + - timezonefinder=6.5.2 - scikit-image=0.24.0 - exactextract=0.2.0 - pip: diff --git a/notebooks/layers/era5.ipynb b/notebooks/layers/era5.ipynb new file mode 100644 index 0000000..05b3b28 --- /dev/null +++ b/notebooks/layers/era5.ipynb @@ -0,0 +1,124 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Setup" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "sys.dont_write_bytecode=True\n", + "\n", + "%load_ext autoreload\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import geopandas as gpd" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# # update the wd path to be able to laod the module\n", + "os.chdir('../..')\n", + "os.getcwd()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Get Area of Interest" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# load boundary from s3\n", + "boundary_path = 'https://cities-indicators.s3.eu-west-3.amazonaws.com/data/boundaries/boundary-BRA-Salvador-ADM4union.geojson'\n", + "city_gdf = gpd.read_file(boundary_path, driver='GeoJSON')\n", + "city_gdf.head()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get area in sqare km\n", + "city_gdf.to_crs(epsg=3857).area / 10**6" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get Layer" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%autoreload\n", + "from city_metrix.layers import Era5HottestDay\n", + "\n", + "hottest_day = Era5HottestDay().get_data(city_gdf.total_bounds)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "hottest_day" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "cities-cif", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/notebooks/layers/era_5_hottest_day.ipynb b/notebooks/layers/era_5_hottest_day.ipynb new file mode 100644 index 0000000..fee05e5 --- /dev/null +++ b/notebooks/layers/era_5_hottest_day.ipynb @@ -0,0 +1,385 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Setup" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "import geopandas as gpd" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "'/home/weiqi_tori/GitHub/wri/cities-cif'" + ] + }, + "execution_count": 2, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# # update the wd path to be able to laod the module\n", + "os.chdir('../..')\n", + "os.getcwd()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Get Area of Interest" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/weiqi_tori/anaconda3/envs/fenv/lib/python3.10/site-packages/pyogrio/raw.py:196: RuntimeWarning: driver GeoJSON does not support open option DRIVER\n", + " return ogr_read(\n" + ] + }, + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
geo_idgeo_levelgeo_namegeo_parent_namecreation_dategeometry
0BRA-Salvador_ADM4-union_1ADM4-unionBRA-SalvadorBRA-Salvador2022-08-03MULTIPOLYGON (((-38.50135 -13.01134, -38.5014 ...
\n", + "
" + ], + "text/plain": [ + " geo_id geo_level geo_name geo_parent_name \\\n", + "0 BRA-Salvador_ADM4-union_1 ADM4-union BRA-Salvador BRA-Salvador \n", + "\n", + " creation_date geometry \n", + "0 2022-08-03 MULTIPOLYGON (((-38.50135 -13.01134, -38.5014 ... " + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# load boundary from s3\n", + "boundary_path = 'https://cities-indicators.s3.eu-west-3.amazonaws.com/data/boundaries/boundary-BRA-Salvador-ADM4union.geojson'\n", + "city_gdf = gpd.read_file(boundary_path, driver='GeoJSON')\n", + "city_gdf.head()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Get Layer" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Authenticating to GEE with configured credentials file.\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "2024-09-20 17:16:16,081 INFO Request ID is 2b411663-6790-42c7-ab4d-f3ae38c99678\n", + "2024-09-20 17:16:16,431 INFO status has been updated to accepted\n", + "2024-09-20 17:16:38,041 INFO status has been updated to successful\n" + ] + }, + { + "data": { + "application/vnd.jupyter.widget-view+json": { + "model_id": "550d881c31b540dcbb82290f7e9323c8", + "version_major": 2, + "version_minor": 0 + }, + "text/plain": [ + "f4e4b5329c3eea8496d48601d058dbae.nc: 0%| | 0.00/100k [00:00\n", + " .geemap-dark {\n", + " --jp-widgets-color: white;\n", + " --jp-widgets-label-color: white;\n", + " --jp-ui-font-color1: white;\n", + " --jp-layout-color2: #454545;\n", + " background-color: #383838;\n", + " }\n", + "\n", + " .geemap-dark .jupyter-button {\n", + " --jp-layout-color3: #383838;\n", + " }\n", + "\n", + " .geemap-colab {\n", + " background-color: var(--colab-primary-surface-color, white);\n", + " }\n", + "\n", + " .geemap-colab .jupyter-button {\n", + " --jp-layout-color3: var(--colab-primary-surface-color, white);\n", + " }\n", + " \n", + " " + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAHFCAYAAAAExnZzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/TGe4hAAAACXBIWXMAAA9hAAAPYQGoP6dpAABxrElEQVR4nO3deVxUVf8H8M+AMqwzisgmCOaOiqmUQhaIgWKP4pK5lGIauZEhmaaWoAWaT25PlpUZamlYKWalJKZQpiaoGKKhEggZiKgsigLC+f1h3J8j2wwzLE6f9+91Xz/vueee5Yrxfc4591yZEEKAiIiIiOpk0NQNICIiInpYMHAiIiIiUhMDJyIiIiI1MXAiIiIiUhMDJyIiIiI1MXAiIiIiUhMDJyIiIiI1MXAiIiIiUhMDJyIiIiI1MXCif4W4uDjIZDJ88803Td0UnYqKisKjjz4KY2Nj2NvbIzg4GDdv3qzzvs2bN0Mmk0lHXl6edC0lJQWzZs2Cu7s7zMzMIJPJEBcXV205rVq1ksoICgrSVbcajJeXF7y8vOp1r7OzM/7zn//UmS8jIwMymQybN2+uVz0N7cMPP6y2befPn8e8efPQr18/tGrVCpaWlnjiiSf07t8MkbYYOBE9pLZt24YJEybgsccew759+xAaGorNmzdj9OjRapexa9cuHD16FK1atZLSEhMTsXv3blhaWmLw4MG13n/gwAEcPXq0vl1odB9++CE+/PDDpm5Gk6opcNq/fz9++OEHjBkzBl9//TW2bduGzp07Y+zYsVi2bFnjN5SomWrR1A0g0ie3b9+GsbExZDJZg9ZTXl6O119/Hb6+vti4cSMAYNCgQbCwsMDzzz+Pffv2wc/Pr85y+vTpA2dnZ5W0SZMmISAgAADwzTff4Lvvvqvxfjc3t/p3ohEVFxfD1NQULi4uTd2UZmv8+PGYPXu2ys+un58f8vLy8O6772LBggWQy+VN2EKi5oEjTqRzYWFhkMlkSElJwYQJE6BUKmFjY4OpU6eioKBAylfblIZMJkNYWFiVMn///XeMHTsWSqUSlpaWCAkJwd27d5GamoqhQ4fCwsICzs7OWLlyZbVtu3PnDkJCQmBrawsTExN4enri1KlTVfIlJiZixIgRsLS0hLGxMfr06YOvvvpKJU/ldNf+/fsxdepUtG3bFqampigpKanfg9PAsWPHkJ2djRdffFElfezYsTA3N0d0dHS9yzYwaLz/LIwcORJOTk6oqKiocq1///7o27evdP7BBx/gqaeegrW1NczMzNCrVy+sXLkSZWVlKvd5eXmhZ8+e+Pnnn+Hh4QFTU1NMnTpVuvbgVN3SpUvRv39/WFpaQqFQoG/fvti0aRNq+v55dHQ0XF1dYWxsjEceeQT/+9//1OrrhQsXMHHiRFhbW0Mul6N79+744IMP1LpXV5ydnZGSkoL4+HhpirUycLaysqo24H/88cdRXFyM69evS2lTpkyBubk5/vjjDwwZMgRmZmaws7PDihUrANz7+Rw4cCDMzMzQpUsXbNmypVH6R9QYOOJEDWbMmDEYN24cpk2bhuTkZCxcuBAA8Nlnn9W7zOeeew4vvPACpk+fjtjYWOkX54EDBzBr1izMmzcP27dvx4IFC9CpU6cq01aLFi1C37598emnn6KgoABhYWHw8vLCqVOn8MgjjwAADh06hKFDh6J///746KOPoFQqERUVhXHjxqG4uBhTpkxRKXPq1Kl45pln8Pnnn+PWrVto2bJlje2/e/euWv00NDSsddTqzJkzAABXV1eV9JYtW6Jbt27S9eZu6tSp8Pf3x8GDB/H0009L6X/88QeOHz+uEpSkpaVh4sSJ6NChA4yMjHD69GmEh4fjjz/+qPIzlZ2djRdeeAHz589HRERErcFgRkYGpk+fjvbt2wO490v/lVdeweXLl7FkyRKVvElJSQgODkZYWBhsbW2xbds2vPrqqygtLcW8efNqrOPs2bPw8PBA+/btsWrVKtja2uLHH3/EnDlzkJeXh9DQ0FqfU3l5eY2B3P0MDAxq7Wt0dDSeffZZKJVKacqyrlGkQ4cOoW3btrC2tlZJLysrw+jRozFjxgy8/vrr2L59OxYuXIjCwkLs3LkTCxYsgIODA95//31MmTIFPXv2RL9+/ersA1GzJ4h0LDQ0VAAQK1euVEmfNWuWMDY2FhUVFUIIIdLT0wUAERkZWaUMACI0NLRKmatWrVLJ9+ijjwoAYteuXVJaWVmZaNu2rRg9erSUdujQIQFA9O3bV6pfCCEyMjJEy5YtxUsvvSSldevWTfTp00eUlZWp1PWf//xH2NnZifLyciGEEJGRkQKAmDx5sppP5l6/1Dmqeyb3Cw8PFwBEdnZ2lWu+vr6iS5cutd5f2fb09PRa83399dcCgDh06FCd/Zo9e3ateapTVlYmbGxsxMSJE1XS58+fL4yMjEReXl6195WXl4uysjKxdetWYWhoKK5fvy5d8/T0FADETz/9VOU+T09P4enpWWN7KstdtmyZaNOmjcrPipOTk5DJZCIpKUnlHh8fH6FQKMStW7eEENX/XA8ZMkQ4ODiIgoIClXuDgoKEsbGxSvurU9mnuo6AgIBayxFCiB49etT6DO63ceNGAUCsW7dOJT0gIEAAEDt37pTSKv/dARAnT56U0q9duyYMDQ1FSEiIWnUSNXcccaIGM2LECJVzV1dX3LlzB7m5ubCxsalXmQ++1dS9e3ecPn1aZT1PixYt0KlTJ1y6dKnK/RMnTlQZyXFycoKHhwcOHToEALh48SL++OMPvPfeewBUR4iGDRuG77//HqmpqejevbuUPmbMGLXbn5CQoFa+Dh06qJWvplGphl5jpSstWrTACy+8gA8++AAFBQVQKpUoLy/H559/Dn9/f7Rp00bKe+rUKYSGhuLXX39VmTYC7r0R1r9/f+m8devW8Pb2VqsNBw8eREREBBISElBYWKhy7cGf1R49eqB3794qeSZOnIjY2FicPHkSAwcOrFL+nTt38NNPP2HmzJkwNTWt8jO1fv16HDt2rNY1aR9//DGKiorq7IuVlVWdedS1b98+zJ49G88++yxeeeWVKtdlMhmGDRsmnVf+u2vRogX69OkjpVtaWsLa2rraf49EDyMGTtRg7v+lB/z/lMDt27frXaalpaXKuZGREUxNTWFsbFwl/cFfggBga2tbbdrp06cBAFeuXAEAzJs3r8apl/tf3QcAOzs7tdv/6KOPqpXP0NCw1uuVz/batWtVgtDr169XeU7N2dSpU7Fq1SpERUVh+vTp+PHHH6us38rMzMSTTz6Jrl27Yt26dXB2doaxsTGOHz+O2bNnV/mZUvfv5Pjx4/D19YWXlxc2btwIBwcHGBkZYffu3QgPD69Sbk0/P8C9v4vqXLt2DXfv3sX777+P999/v9o8D/5MPahTp05qT9Xpwo8//ojRo0fDx8cH27ZtqzYQr+nfXXU/e0ZGRrhz545O2kbU1Bg4UZOp/I/ug4upa/oFpAs5OTnVplUGIpX/i33hwoU1vtbftWtXlXNNRndqW/90v8jIyCprqe7Xq1cvAEBycrLKm2J3797FH3/8gQkTJqjdpqbm4uKCxx9/HJGRkZg+fToiIyNhb28PX19fKc/u3btx69Yt7Nq1C05OTlJ6UlJStWWq+3cSFRWFli1b4vvvv1cJAnbv3l1t/pp+foCq/0OhUuvWrWFoaIhJkyZh9uzZ1eapa4Rx8ODBiI+PrzUPAAQEBGi9f9SPP/6IkSNHwtPTEzt37oSRkZFW5RHpGwZO1GRsbGxgbGyM33//XSX922+/bbA6v/zyS4SEhEi/WC9duoQjR45g8uTJAO4FRZ07d8bp06cRERGh8/p1NVXXv39/2NnZYfPmzRg3bpyU/s033+DmzZsa7eXUHLz44ouYOXMmDh8+jO+++w4hISEqo26Vf1/3L2QWQkhbMdSXTCZDixYtVOq6ffs2Pv/882rzp6Sk4PTp0yrTddu3b4eFhYXKG4D3MzU1xaBBg3Dq1Cm4urrWKxDR5VSdXC6vcdR3//79GDlyJAYOHIjdu3dz+wGiajBwoiYjk8nwwgsv4LPPPkPHjh3Ru3dvHD9+HNu3b2+wOnNzczFq1CgEBgaioKAAoaGhMDY2lt74A+79kvLz88OQIUMwZcoUtGvXDtevX8e5c+dw8uRJfP311/WuX1f7HhkaGmLlypWYNGkSpk+fjgkTJuDChQuYP38+fHx8MHTo0HqXXVxcjL179wK494YZAMTHxyMvLw9mZmZq7Q8F3Pv79fT0rHHX8ftNmDABISEhmDBhAkpKSqqMtvn4+MDIyAgTJkzA/PnzcefOHWzYsAE3btzQqG8PeuaZZ7B69WpMnDgRL7/8Mq5du4b33nuvxoDB3t4eI0aMQFhYGOzs7PDFF18gNjYW7777LkxNTWusZ926dRg4cCCefPJJzJw5E87OzigqKsLFixfx3Xff4eDBg7W288FRTm306tULUVFR2LFjBx555BEYGxujV69eOHz4MEaOHAlbW1ssWrSoymiei4sLFAqFztpB9LBi4ERNatWqVQCAlStX4ubNm/D29sb3339fZVNGXalcBPziiy+isLAQjz/+OKKiotCxY0cpz6BBg3D8+HGEh4cjODgYN27cQJs2beDi4oLnnnuuQdpVHy+88AIMDQ2xYsUKbN68GZaWlpg8eTLCw8O1Kjc3Nxdjx45VSavcU8vJyQkZGRl1llH52Rd11xoplUqMGjUK27dvxxNPPIEuXbqoXO/WrRt27tyJN998E6NHj0abNm0wceJEhISEqB3IVcfb2xufffYZ3n33XQwfPhzt2rVDYGAgrK2tMW3atCr5H330Ubz44osIDQ3FhQsXYG9vj9WrV2Pu3Lm11uPi4oKTJ0/i7bffxptvvonc3Fy0atUKnTt3Vllg3RiWLl2K7OxsBAYGoqioSPo7PXDgAG7fvo2MjIxqF9YfOnSo3p+rIdInMqHOikMi0iubN2/Giy++iIsXL8LJyQktWtTvf0NV7i/UsmVLzJ49G+vXrwcA7N27F//5z39w+vRpaT0WEZE+4M7hRP9inTp1QsuWLet8q6smbdq0qXbB+6FDhzB+/HgGTUSkdzjiRPQvdO3aNaSnp0vnjz76aL1GnZKSkqR9iaytraXdt4mI9BUDJyIiIiI1caqOiIiISE0MnIiIiIjUxMCJiIiISE3cx0kHKioq8Pfff8PCwuKh+bgqERE1DSEEioqKYG9vr7PvC1bnzp07KC0t1bocIyOjKt8l/Ddj4KQDf//9NxwdHZu6GURE9BDJysqCg4NDg5R9584dtDExRTG0f//L1tYW6enpDJ7+wcBJBywsLAAA31vZwawB/9cD0cNma25BUzeBqNkphcA23JJ+dzRIHaWlKIbAJJjBCPWfCSmFwOc5OSgtLWXg9A8GTjpQOT1nZmAAcwZORBJt/oNNpO8aY2mHEWT8d6hjDJyIiIj0lAFkMNAiQDPgTo9VMHAiIiLSUwbQ7vV5zqFUxWdCREREpCaOOBEREekpmQww0GKJkwyADl7M0ysMnIiIiPQUp+p0j8+EiIiISE0MnIiIiPSUgUym9aGJDRs2wNXVFQqFAgqFAu7u7ti3b590XQiBsLAw2Nvbw8TEBF5eXkhJSalSztGjR+Ht7Q0zMzO0atUKXl5euH37do31Ll++HI899hgsLCxgbW2NkSNHIjU1VSXPlClTIJPJVI4BAwZo1D+AgRMREZHeMtDBoQkHBwesWLECiYmJSExMhLe3N/z9/aXgaOXKlVi9ejXWr1+PhIQE2NrawsfHB0VFRVIZR48exdChQ+Hr64vjx48jISEBQUFBtX6eJj4+HrNnz8axY8cQGxuLu3fvwtfXF7du3VLJN3ToUGRnZ0vH3r17NewhIBNCcNmXlgoLC6FUKnHIuh03wCS6z6c5+U3dBKJmpxQCkbiJgoICKBSKBqmj8vfSHAMF5Frs41QiBP5XUahVWy0tLfHf//4XU6dOhb29PYKDg7FgwYJ75ZeUwMbGBu+++y6mT58OABgwYAB8fHzw9ttv17vdV69ehbW1NeLj4/HUU08BuDfilJ+fj927d9e7XIAjTkRERHrLQKb9UV/l5eWIiorCrVu34O7ujvT0dOTk5MDX11fKI5fL4enpiSNHjgAAcnNz8dtvv8Ha2hoeHh6wsbGBp6cnDh8+rFHdBQX3PvdkaWmpkh4XFwdra2t06dIFgYGByM3N1bhfDJyIiIj0lK6m6goLC1WOkpKSGutMTk6Gubk55HI5ZsyYgejoaLi4uCAnJwcAYGNjo5LfxsZGuvbnn38CAMLCwhAYGIiYmBj07dsXgwcPxoULF9TqsxACISEhGDhwIHr27Cml+/n5Ydu2bTh48CBWrVqFhIQEeHt719qX6nA7AiIiIj1VuQi63vf/8/8dHR1V0kNDQxEWFlbtPV27dkVSUhLy8/Oxc+dOBAQEID4+XqVN9xNCSGkVFRUAgOnTp+PFF18EAPTp0wc//fQTPvvsMyxfvrzONgcFBeH333+vMko1btw46c89e/aEm5sbnJyc8MMPP2D06NF1lluJgRMRERHVKisrS2WNk1wurzGvkZEROnXqBABwc3NDQkIC1q1bJ61rysnJgZ2dnZQ/NzdXGoWqTHdxcVEps3v37sjMzKyzna+88gr27NmDn3/+GQ4ODrXmtbOzg5OTk9ojWZU4VUdERKSndDVVV7m9QOVRW+D0ICEESkpK0KFDB9ja2iI2Nla6Vlpaivj4eHh4eAAAnJ2dYW9vX2UrgfPnz8PJyanWOoKCgrBr1y4cPHgQHTp0qLNd165dQ1ZWlkoQpw6OOBEREekpbRd4azq6smjRIvj5+cHR0RFFRUWIiopCXFwcYmJiIJPJEBwcjIiICHTu3BmdO3dGREQETE1NMXHiRAD3pvFef/11hIaGonfv3nj00UexZcsW/PHHH/jmm2+kegYPHoxRo0YhKCgIADB79mxs374d3377LSwsLKQ1U0qlEiYmJrh58ybCwsIwZswY2NnZISMjA4sWLYKVlRVGjRqlUR8ZOBEREZFOXLlyBZMmTUJ2djaUSiVcXV0RExMDHx8fAMD8+fNx+/ZtzJo1Czdu3ED//v2xf/9+WFhYSGUEBwfjzp07mDt3Lq5fv47evXsjNjYWHTt2lPKkpaUhLy9POt+wYQMAwMvLS6U9kZGRmDJlCgwNDZGcnIytW7ciPz8fdnZ2GDRoEHbs2KFStzq4j5MOcB8noupxHyeiqhpzH6c3WiphrMXi8DtCYEVZQYO29WHDESciIiI9VZ/Ppqjcr8O26As+EyIiIiI1ccSJiIhIT9Xne3MP3k+qGDgRERHpqcZ+q+7fgM+EiIiISE0ccSIiItJTnKrTPQZOREREesoAMhiAb9XpEgMnIiIiPcU1TrrHZ0JERESkJo44ERER6SmucdI9Bk5ERER6ilN1usdnQkRERKQmjjgRERHpKRmg1Vt1MgjdNUZPMHAiIiLSU5yq0z0+EyIiIiI1ccSJiIhIT/GtOt1j4ERERKSnOFWne3wmRERERGriiBMREZGe0v5bdVoMV+mph2bEKTw8HB4eHjA1NUWrVq2qXD99+jQmTJgAR0dHmJiYoHv37li3bl2tZWZkZEAmk1V7fP311w3UEyIiosZROVWnzUGqHpoRp9LSUowdOxbu7u7YtGlTlesnTpxA27Zt8cUXX8DR0RFHjhzByy+/DENDQwQFBVVbpqOjI7Kzs1XSPvnkE6xcuRJ+fn4N0g8iIqLGIvvn0OZ+UvXQBE5Lly4FAGzevLna61OnTlU5f+SRR3D06FHs2rWrxsDJ0NAQtra2KmnR0dEYN24czM3NtW80ERER6ZWHJnCqj4KCAlhaWqqd/8SJE0hKSsIHH3xQa76SkhKUlJRI54WFhfVuIxERUUPhW3W6p7fP5OjRo/jqq68wffp0te/ZtGkTunfvDg8Pj1rzLV++HEqlUjocHR21bS4REZHOVS4O1+YgVU0aOIWFhdW4OLvySExM1LjclJQU+Pv7Y8mSJfDx8VHrntu3b2P79u2YNm1anXkXLlyIgoIC6cjKytK4jURERPTwadKpuqCgIIwfP77WPM7OzhqVefbsWXh7eyMwMBBvvvmm2vd98803KC4uxuTJk+vMK5fLIZfLNWoXERFRY+NUne41aeBkZWUFKysrnZWXkpICb29vBAQEIDw8XKN7N23ahBEjRqBt27Y6aw8REVFTkkG74IcTdVU9NMFkZmYmkpKSkJmZifLyciQlJSEpKQk3b94EcC9oGjRoEHx8fBASEoKcnBzk5OTg6tWrUhmXL19Gt27dcPz4cZWyL168iJ9//hkvvfRSo/aJiIiIHi4PzVt1S5YswZYtW6TzPn36AAAOHToELy8vfP3117h69Sq2bduGbdu2SfmcnJyQkZEBACgrK0NqaiqKi4tVyv7ss8/Qrl07+Pr6NnxHiIiIGgn3cdI9mRBCNHUjHnaFhYVQKpU4ZN0O5gYPzSAeUYP7NCe/qZtA1OyUQiASN1FQUACFQtEgdVT+XvpMaQVTWf1/LxWLCkwtyGvQtj5s+FueiIiISE0MnIiIiPSUTAeHJjZs2ABXV1coFAooFAq4u7tj37590nUhBMLCwmBvbw8TExN4eXkhJSWlSjlHjx6Ft7c3zMzM0KpVK3h5eeH27du11v3hhx+iQ4cOMDY2Rr9+/fDLL7+oXFe37rowcCIiItJTjR04OTg4YMWKFUhMTERiYiK8vb3h7+8vBSgrV67E6tWrsX79eiQkJMDW1hY+Pj4oKiqSyjh69CiGDh0KX19fHD9+HAkJCQgKCoJBLUthduzYgeDgYCxevBinTp3Ck08+CT8/P2RmZkp51KlbHVzjpANc40RUPa5xIqqqMdc4bdHBGqcALdc4WVpa4r///S+mTp0Ke3t7BAcHY8GCBQDufcLMxsYG7777rvSljwEDBsDHxwdvv/222nX0798fffv2xYYNG6S07t27Y+TIkVi+fDmEEGrVrQ7+liciIqJaFRYWqhz3f6+1JuXl5YiKisKtW7fg7u6O9PR05OTkqLzBLpfL4enpiSNHjgAAcnNz8dtvv8Ha2hoeHh6wsbGBp6cnDh8+XGM9paWlOHHiRJU34319faVy1albXQyciIiI9FRdnzVT5wAAR0dHlW+0Ll++vMY6k5OTYW5uDrlcjhkzZiA6OhouLi7IyckBANjY2Kjkt7Gxka79+eefAO59ki0wMBAxMTHo27cvBg8ejAsXLlRbX15eHsrLy2stV5261fXQ7ONEREREmtHVPk5ZWVkqU3W1fXasa9euSEpKQn5+Pnbu3ImAgADEx8f/f5ky1RYJIaS0iooKAMD06dPx4osvAri3b+NPP/2Ezz77rNaArbZyNclTFwZOREREVKvKt+TUYWRkhE6dOgEA3NzckJCQgHXr1klri3JycmBnZyflz83NlUaCKtNdXFxUyuzevbvKQu/7WVlZwdDQsMrI0f3l2tra1lm3ujhVR0REpKcMdHBoSwiBkpISdOjQAba2toiNjZWulZaWIj4+Hh4eHgAAZ2dn2NvbIzU1VaWM8+fPw8nJqdryjYyM0K9fP5VyASA2NlYqV5261cURJyIiIj0lk9076n2/hvkXLVoEPz8/ODo6oqioCFFRUYiLi0NMTAxkMhmCg4MRERGBzp07o3PnzoiIiICpqSkmTpz4T3tleP311xEaGorevXvj0UcfxZYtW/DHH3/gm2++keoZPHgwRo0ahaCgIABASEgIJk2aBDc3N7i7u+OTTz5BZmYmZsyYIZVbV93qYuBEREREOnHlyhVMmjQJ2dnZUCqVcHV1RUxMDHx8fAAA8+fPx+3btzFr1izcuHED/fv3x/79+2FhYSGVERwcjDt37mDu3Lm4fv06evfujdjYWHTs2FHKk5aWhry8POl83LhxuHbtGpYtW4bs7Gz07NkTe/fuVRmlUqdudXAfJx3gPk5E1eM+TkRVNeY+Tl+2ttZ6H6cJN3L5rbr7cMSJiIhIT+nqrTr6fxweISIiIlITR5yIiIj0FEecdI+BExERkZ4yAGCgRfRjwFXQVTBwIiIi0lOyf/5Pm/tJFdc4EREREamJI05ERER6jGNGusXAiYiISE9pvXM4o64qOFVHREREpCaOOBEREekpbkegewyciIiI9JQBZDDQIvzR5l59xak6IiIiIjVxxImIiEhPcapO9xg4ERER6Sm+Vad7nKojIiIiUhNHnIiIiPQUp+p0j4ETERGRnuK36nSPgRMREZGeMpDdO7S5n1RxjRMRERGRmjjiREREpKe4xkn3GDgRERHpKQZOusepOiIiIiI1ccSJiIhIT/GtOt1j4ERERKSnuHO47nGqjoiIiEhNHHEiIiLSUwbQboSEoytVMXAiIiLSU3yrTvcYOBEREekrmQwyLnLSKY7CEREREamJI05ERER6ilN1uscRJyIiIj0l08GhiQ0bNsDV1RUKhQIKhQLu7u7Yt2+fdF0IgbCwMNjb28PExAReXl5ISUlRKcPLywuyf6YYK4/x48fXWq+zs3OVe2QyGWbPni3lmTJlSpXrAwYM0LCHDJyIiIhIRxwcHLBixQokJiYiMTER3t7e8Pf3l4KjlStXYvXq1Vi/fj0SEhJga2sLHx8fFBUVqZQTGBiI7Oxs6fj4449rrTchIUElf2xsLABg7NixKvmGDh2qkm/v3r0a95FTdURERHpKpuXicE3vHT58uMp5eHg4NmzYgGPHjsHFxQVr167F4sWLMXr0aADAli1bYGNjg+3bt2P69OnSfaamprC1tVW73rZt26qcr1ixAh07doSnp6dKulwu16jc6nDEiYiISE8ZyLQ/AKCwsFDlKCkpqbPu8vJyREVF4datW3B3d0d6ejpycnLg6+sr5ZHL5fD09MSRI0dU7t22bRusrKzQo0cPzJs3r8qIVG1KS0vxxRdfYOrUqVUCv7i4OFhbW6NLly4IDAxEbm6u2uVW4ogTERER1crR0VHlPDQ0FGFhYdXmTU5Ohru7O+7cuQNzc3NER0fDxcVFCo5sbGxU8tvY2ODSpUvS+fPPP48OHTrA1tYWZ86cwcKFC3H69Glp+q0uu3fvRn5+PqZMmaKS7ufnh7Fjx8LJyQnp6el466234O3tjRMnTkAul6tVNsDAiYiISG/JDGSQGWgxVffP8vCsrCwoFAopvbZAo2vXrkhKSkJ+fj527tyJgIAAxMfH/3+ZD4wCCSFU0gIDA6U/9+zZE507d4abmxtOnjyJvn371tnmTZs2wc/PD/b29irp48aNUynXzc0NTk5O+OGHH6SpQ3UwcCIiItJTuvrIb+VbcuowMjJCp06dAABubm5ISEjAunXrsGDBAgBATk4O7OzspPy5ublVRqHu17dvX7Rs2RIXLlyoM3C6dOkSDhw4gF27dtXZTjs7Ozg5OeHChQvqdEvCNU5ERETUYIQQKCkpkabf7p9yKy0tRXx8PDw8PGq8PyUlBWVlZSrBVk0iIyNhbW2NZ555ps68165dQ1ZWllrl3o8jTkRERHpKVyNO6lq0aBH8/Pzg6OiIoqIiREVFIS4uDjExMZDJZAgODkZERAQ6d+6Mzp07IyIiAqamppg4cSIAIC0tDdu2bcOwYcNgZWWFs2fP4rXXXkOfPn3wxBNPSPUMHjwYo0aNQlBQkJRWUVGByMhIBAQEoEUL1fDm5s2bCAsLw5gxY2BnZ4eMjAwsWrQIVlZWGDVqlEZ9ZOBERESkpxp7O4IrV65g0qRJyM7OhlKphKurK2JiYuDj4wMAmD9/Pm7fvo1Zs2bhxo0b6N+/P/bv3w8LCwsA96b5fvrpJ6xbtw43b96Eo6MjnnnmGYSGhsLQ0FCqJy0tDXl5eSp1HzhwAJmZmZg6dWqVdhkaGiI5ORlbt25Ffn4+7OzsMGjQIOzYsUOqW+1nIoQQGt1BVRQWFkKpVOKQdTuYG3D2k6jSpzn5Td0EomanFAKRuImCggK11w1pqvL30tH27bX6vXSzogLumZkN2taHjVojTuqsYr+fTCbDnj170K5du3o1ioiIiKg5UitwSkpKwmuvvQZzc/M68wohsGLFCrU2xyIiIqKG09hTdf8Gaq9xev3112Ftba1W3lWrVtW7QURERKQbjb04/N9ArcApPT29yndganP27NkqG08RERERPezUCpycnJw0KvTBrdmJiIio8RnIZDDQYthIm3v1Vb22I8jPz8fx48eRm5uLiooKlWuTJ0/WScOIiIhIO5yq0z2NA6fvvvsOzz//PG7dugULCwuVhWMymYyBExEREektjTd3eO211zB16lQUFRUhPz8fN27ckI7r1683RBuJiIioHmSQSW/W1esAh5wepPGI0+XLlzFnzhyYmpo2RHuIiIhIR2QG9456388tsqvQ+HEOGTIEiYmJDdEWIiIiomZNrRGnPXv2SH9+5pln8Prrr+Ps2bPo1asXWrZsqZJ3xIgRum3hP8LDw/HDDz8gKSkJRkZGyM/PV7l++vRprFixAocPH0ZeXh6cnZ0xY8YMvPrqq7WWm5OTg9dffx2xsbEoKipC165dsWjRIjz77LMN0g8iIqJGo+UGmFwdXpVagdPIkSOrpC1btqxKmkwmQ3l5udaNqk5paSnGjh0Ld3d3bNq0qcr1EydOoG3btvjiiy/g6OiII0eO4OWXX4ahoaHK15MfNGnSJBQUFGDPnj2wsrLC9u3bMW7cOCQmJqJPnz4N0hciIqLGwLfqdE+twOnBLQeawtKlSwEAmzdvrvb6g19DfuSRR3D06FHs2rWr1sDp6NGj2LBhAx5//HEAwJtvvok1a9bg5MmTDJyIiOihdi9w0uaTKzpsjJ7QeI3T1q1bq/0OXWlpKbZu3aqTRulKQUEBLC0ta80zcOBA7NixA9evX0dFRQWioqJQUlICLy+vGu8pKSlBYWGhykFERET6T+PA6cUXX0RBQUGV9KKiIrz44os6aZQuHD16FF999RWmT59ea74dO3bg7t27aNOmDeRyOaZPn47o6Gh07NixxnuWL18OpVIpHdwpnYiImqPKqTptDlKlceAkhKh22O+vv/6CUqnUqKywsLA695Cozxt8KSkp8Pf3x5IlS+Dj41Nr3jfffBM3btzAgQMHkJiYiJCQEIwdOxbJyck13rNw4UIUFBRIR1ZWlsZtJCIiamiVn1zR5iBVau/j1KdPHymYGTx4MFq0+P9by8vLkZ6ejqFDh2pUeVBQEMaPH19rHmdnZ43KPHv2LLy9vREYGIg333yz1rxpaWlYv349zpw5gx49egAAevfujV9++QUffPABPvroo2rvk8vlkMvlGrWLiIiIHn5qB06Vb9YlJSVhyJAhMDc3l64ZGRnB2dkZY8aM0ahyKysrWFlZaXRPbVJSUuDt7Y2AgACEh4fXmb+4uBgAYGCgOvBmaGjYLBbEExERaYNv1eme2oFTaGgoysvL4eTkhCFDhsDOzq4h21VFZmYmrl+/jszMTJSXlyMpKQkA0KlTJ5ibmyMlJQWDBg2Cr68vQkJCkJOTA+BeENS2bVsA93Y9Hzx4MLZu3YrHH38c3bp1Q6dOnTB9+nS89957aNOmDXbv3o3Y2Fh8//33jdo/IiIiXaucKdLmflKl0RonQ0NDzJgxA3fu3Gmo9tRoyZIl6NOnD0JDQ3Hz5k306dMHffr0kdZAff3117h69Sq2bdsGOzs76XjsscekMsrKypCamiqNNLVs2RJ79+5F27ZtMXz4cLi6umLr1q3YsmULhg0b1uh9JCIiouZNJoTQ6Es0jz32GFasWIHBgwc3VJseOoWFhVAqlThk3Q7mBlp8FIhIz3yak9/UTSBqdkohEImbKCgogEKhaJA6Kn8vne3ZCRaGhvUup6i8HC5nLjZoWx82Gv+WDw8Px7x58/D9998jOzub+xkRERE1U9yOQPfUXuNUqfLNuREjRqjMfVZuU9BQn1whIiIiamoaB06HDh1qiHYQERGRjskMZJAZaLE4XHDI6UEaB06enp4N0Q4iIiLSMW5HoHsaB04AkJ+fj02bNuHcuXOQyWRwcXHB1KlTNd45nIiIiBqOtrt/c+fwqjReHJ6YmIiOHTtizZo1uH79OvLy8rB69Wp07NgRJ0+ebIg2EhERETULGo84zZ07FyNGjMDGjRulz67cvXsXL730EoKDg/Hzzz/rvJFERESkOU7V6Z7GgVNiYqJK0AQALVq0wPz58+Hm5qbTxhEREVH9cedw3dN4qk6hUCAzM7NKelZWFiwsLHTSKCIiIqLmSOPAady4cZg2bRp27NiBrKws/PXXX4iKisJLL72ECRMmNEQbiYiIqB5k0HIDTA3r27BhA1xdXaFQKKBQKODu7o59+/ZJ14UQCAsLg729PUxMTODl5YWUlBSVMry8vKSRsspj/PjxtdYbFhZW5R5bW1uVPOrUrQ6Np+ree+89yGQyTJ48GXfv3gVw75tvM2fOxIoVKzRuABERETWMxp6qc3BwwIoVK9CpUycAwJYtW+Dv749Tp06hR48eWLlyJVavXo3NmzejS5cueOedd+Dj44PU1FSVWavAwEAsW7ZMOjcxMamz7h49euDAgQPSueEDn5pRt+66aBw4GRkZYd26dVi+fDnS0tIghECnTp1gamqqaVFERESkR4YPH65yHh4ejg0bNuDYsWNwcXHB2rVrsXjxYowePRrAvcDKxsYG27dvx/Tp06X7TE1Nq4wY1aVFixY13iOEULvuutT7i7Smpqbo1asXXF1dGTQRERE1R9p+p+6fAacHv0tbUlJSZ9Xl5eWIiorCrVu34O7ujvT0dOTk5MDX11fKI5fL4enpiSNHjqjcu23bNlhZWaFHjx6YN28eioqK6qzvwoULsLe3R4cOHTB+/Hj8+eef0jVN6q6LxiNOt27dwooVK/DTTz8hNzcXFRUVKtfvbygRERE1HV1N1Tk6Oqqkh4aGIiwsrNp7kpOT4e7ujjt37sDc3BzR0dFwcXGRAhQbGxuV/DY2Nrh06ZJ0/vzzz6NDhw6wtbXFmTNnsHDhQpw+fRqxsbE1trN///7YunUrunTpgitXruCdd96Bh4cHUlJS0KZNG+Tk5KhVtzo0DpxeeuklxMfHY9KkSbCzs+OrikRERHouKysLCoVCOpfL5TXm7dq1K5KSkpCfn4+dO3ciICAA8fHx0vUH4wYhhEpaYGCg9OeePXuic+fOcHNzw8mTJ9G3b99q6/Tz85P+3KtXL7i7u6Njx47YsmULQkJC1K5bHRoHTvv27cMPP/yAJ554QtNbiYiIqBHJDO4d2twPQHpLTh1GRkbS4nA3NzckJCRg3bp1WLBgAQAgJycHdnZ2Uv7c3NwqI0H369u3L1q2bIkLFy7UGDg9yMzMDL169cKFCxcAQFr7pGnd1dH4cbZu3RqWlpaa3kZERESN7MFX9OtzaEsIgZKSEmn67f4pt9LSUsTHx8PDw6PG+1NSUlBWVqYS8NSlpKQE586dk+6pb93V0Thwevvtt7FkyRIUFxdreisRERE1JgOZ9ocGFi1ahF9++QUZGRlITk7G4sWLERcXh+effx4ymQzBwcGIiIhAdHQ0zpw5gylTpsDU1BQTJ04EAKSlpWHZsmVITExERkYG9u7di7Fjx6JPnz4qM12DBw/G+vXrpfN58+YhPj4e6enp+O233/Dss8+isLAQAQEBAKBW3erSeKpu1apVSEtLg42NDZydndGyZUuV6/zQLxER0b/TlStXMGnSJGRnZ0OpVMLV1RUxMTHw8fEBAMyfPx+3b9/GrFmzcOPGDfTv3x/79++X9lEyMjLCTz/9hHXr1uHmzZtwdHTEM888g9DQUJV9mdLS0pCXlyed//XXX5gwYQLy8vLQtm1bDBgwAMeOHYOTk5OUp6661SUTQghNbli6dGmt10NDQzVqgD4oLCyEUqnEIet2MDfQYjKZSM98mpPf1E0ganZKIRCJmygoKFB73ZCmKn8vZQ7sCUULw7pvqKmcu+Vof/hMg7b1YaPxiJO6gdGXX36JESNGwMzMTONGERERkfb4kV/da7DhkenTp+PKlSsNVTwRERFRo9N4xEldGs4AEhERka7VY4F3lftJRYMFTkRERNTEpG+naHE/qeBKZiIiIiI1ccSJiIhIT8kMZJBpMd2mzb36ioETERGRvuJUnc412FSdk5NTlc0xiYiIiB5mGo84ZWVlQSaTwcHBAQBw/PhxbN++HS4uLnj55ZelfGfOnNFdK4mIiEhjMpmWU3UccapC4xGniRMn4tChQwDufWXYx8cHx48fx6JFi7Bs2TKdN5CIiIjqqXKqTpuDVGgcOJ05cwaPP/44AOCrr75Cz549ceTIEWzfvh2bN2/WdfuIiIiovgyg5Ud+m7oDzY/Gj6SsrAxyuRwAcODAAYwYMQIA0K1bN2RnZ+u2dURERETNiMaBU48ePfDRRx/hl19+QWxsLIYOHQoA+Pvvv9GmTRudN5CIiIjqp/JbddocpErjwOndd9/Fxx9/DC8vL0yYMAG9e/cGAOzZs0eawiMiIqJmQKtpOi0/16KnNH6rzsvLC3l5eSgsLETr1q2l9JdffhmmpqY6bRwRERFRc1KvZV9CCJw4cQIff/wxioqKAABGRkYMnIiIiJoTvlWncxqPOF26dAlDhw5FZmYmSkpK4OPjAwsLC6xcuRJ37tzBRx991BDtJCIiIg3JDO4d2txPqjR+JK+++irc3Nxw48YNmJiYSOmjRo3CTz/9pNPGERERETUnGo84HT58GL/++iuMjIxU0p2cnHD58mWdNYyIiIi0xG/V6ZzGgVNFRQXKy8urpP/111+wsLDQSaOIiIhIezIDLT+5wrfqqtB4qs7Hxwdr166VzmUyGW7evInQ0FAMGzZMl20jIiIialY0HnFas2YNBg0aBBcXF9y5cwcTJ07EhQsXYGVlhS+//LIh2khERET1wak6ndM4cLK3t0dSUhK+/PJLnDx5EhUVFZg2bRqef/55lcXiRERE1MS03cSSU3VVaBw4AYCJiQmmTp2KqVOn6ro9REREpCPafjaFn1ypql47NHz++ecYOHAg7O3tcenSJQD3pvC+/fZbnTaOiIiIqDnROHDasGEDQkJC4Ofnhxs3bkhv2LVu3Vpl0TgRERE1MX6rTuc0Dpzef/99bNy4EYsXL0aLFv8/0+fm5obk5GSdNo6IiIi0oe3nVhg4PUjjwCk9PR19+vSpki6Xy3Hr1i2dNIqIiIioOdI4cOrQoQOSkpKqpO/btw8uLi66aBMRERHpQOXicG0OUqXxW3Wvv/46Zs+ejTt37kAIgePHj+PLL7/E8uXL8emnnzZEG4mIiKg+uB2Bzmk84vTiiy8iNDQU8+fPR3FxMSZOnIiPPvoI69atw/jx4xuijURERPQQ2LBhA1xdXaFQKKBQKODu7o59+/ZJ14UQCAsLg729PUxMTODl5YWUlBSVMry8vKqMetUVXyxfvhyPPfYYLCwsYG1tjZEjRyI1NVUlz5QpU6qUO2DAAI37qFHgdPfuXWzZsgXDhw/HpUuXkJubi5ycHGRlZWHatGkaV05EREQNp7Gn6hwcHLBixQokJiYiMTER3t7e8Pf3l4KjlStXYvXq1Vi/fj0SEhJga2sLHx8fFBUVqZQTGBiI7Oxs6fj4449rrTc+Ph6zZ8/GsWPHEBsbi7t378LX17fK2uuhQ4eqlLt3716N+gdoOFXXokULzJw5E+fOnQMAWFlZaVwhERERNZJGnqobPny4ynl4eDg2bNiAY8eOwcXFBWvXrsXixYsxevRoAMCWLVtgY2OD7du3Y/r06dJ9pqamsLW1VbvemJgYlfPIyEhYW1vjxIkTeOqpp6R0uVyuUbnV0Xiqrn///jh16pRWlRIREdHDo7CwUOUoKSmp857y8nJERUXh1q1bcHd3R3p6OnJycuDr6yvlkcvl8PT0xJEjR1Tu3bZtG6ysrNCjRw/MmzevyohUXQoKCgAAlpaWKulxcXGwtrZGly5dEBgYiNzcXI3KBeqxOHzWrFl47bXX8Ndff6Ffv34wMzNTue7q6qpxI4iIiKgB6Ogjv46OjirJoaGhCAsLq/aW5ORkuLu7486dOzA3N0d0dDRcXFyk4MjGxkYlv42NjfQVEgB4/vnn0aFDB9ja2uLMmTNYuHAhTp8+jdjYWLWaLIRASEgIBg4ciJ49e0rpfn5+GDt2LJycnJCeno633noL3t7eOHHiBORyuVplA/UInMaNGwcAmDNnjpQmk8kghIBMJpN2EiciIqKmJTOQQabFVF3lvVlZWVAoFFJ6bYFG165dkZSUhPz8fOzcuRMBAQGIj4///zIfCOQq44dKgYGB0p979uyJzp07w83NDSdPnkTfvn3rbHNQUBB+//13HD58WCW9Mn6pLNfNzQ1OTk744YcfpKlDdWgcOKWnp2t6CxERETUFHY04Vb4lpw4jIyN06tQJwL2viiQkJGDdunVYsGABACAnJwd2dnZS/tzc3CqjUPfr27cvWrZsiQsXLtQZOL3yyivYs2cPfv75Zzg4ONSa187ODk5OTrhw4YJa/aqkceDk5OSk6S1ERET0LyWEQElJiTT9FhsbK32BpLS0FPHx8Xj33XdrvD8lJQVlZWUqwVZ1dbzyyiuIjo5GXFwcOnToUGe7rl27hqysrFrLrY7GgdOePXuqTZfJZDA2NkanTp3UajARERE1MANo+VadZtkXLVoEPz8/ODo6oqioCFFRUYiLi0NMTAxkMhmCg4MRERGBzp07o3PnzoiIiICpqSkmTpwIAEhLS8O2bdswbNgwWFlZ4ezZs3jttdfQp08fPPHEE1I9gwcPxqhRoxAUFAQAmD17NrZv345vv/0WFhYWyMnJAQAolUqYmJjg5s2bCAsLw5gxY2BnZ4eMjAwsWrQIVlZWGDVqlEZ91DhwGjlypLSm6X73r3MaOHAgdu/ejdatW2taPBEREemItp9N0fTeK1euYNKkScjOzoZSqYSrqytiYmLg4+MDAJg/fz5u376NWbNm4caNG+jfvz/2798PCwsLAPem+X766SesW7cON2/ehKOjI5555hmEhobC0NBQqictLQ15eXnS+YYNGwDc2zzzfpGRkZgyZQoMDQ2RnJyMrVu3Ij8/H3Z2dhg0aBB27Ngh1a32MxEPRkB1+Omnn7B48WKEh4fj8ccfBwAcP34cb775Jt566y0olUpMnz4d/fv3x6ZNmzRqzMOqsLAQSqUSh6zbwdxA4x0eiPTWpzn5Td0EomanFAKRuImCggK11w1pqvL30rVpPlAYtax/OaVlaLMptkHb+rDReMTp1VdfxSeffAIPDw8pbfDgwTA2NsbLL7+MlJQUrF27FlOnTtVpQx8GfU/8DIVCs8iVSJ+5tTBq6iYQNTuFhYWItGvfOJXxW3U6p/HwSFpaWrVRp0KhwJ9//gkA6Ny5s8oQmi6Eh4fDw8MDpqamaNWqVZXrp0+fxoQJE+Do6AgTExN0794d69atq7PctLQ0jBo1Cm3btoVCocBzzz2HK1eu6LTtRERETaLyrTptDlKhceDUr18/vP7667h69aqUdvXqVcyfPx+PPfYYAODChQt1vgaoqdLSUowdOxYzZ86s9vqJEyfQtm1bfPHFF0hJScHixYuxcOFCrF+/vsYyb926BV9fX8hkMhw8eBC//vorSktLMXz4cFRUVOi0/URERPTw03iqbtOmTfD394eDgwMcHR0hk8mQmZmJRx55BN9++y0A4ObNm3jrrbd02tClS5cCADZv3lzt9QenBh955BEcPXoUu3btklbdP+jXX39FRkYGTp06JY2iRUZGwtLSEgcPHsTTTz+tuw4QERE1Nh3t40T/T+PAqWvXrjh37hx+/PFHnD9/HkIIdOvWDT4+PjD4Z2H0yJEjdd3OeikoKKjynZr7lZSUQCaTqeyAamxsDAMDAxw+fLjGwKmkpETlOz2FhYW6azQREZHOaDvdxsDpQRoHTsC91xOHDh0KLy8vyOVyrV51bChHjx7FV199hR9++KHGPAMGDICZmRkWLFiAiIgICCGwYMECVFRUIDs7u8b7li9fLo2AERER0b+HxmucKioq8Pbbb6Ndu3YwNzeXPsHy1ltvabz9QFhYmLTHRE1HYmKipk1ESkoK/P39sWTJEmnviOq0bdsWX3/9Nb777juYm5tDqVSioKAAffv2Vdkv4kELFy5EQUGBdGRlZWncRiIiogZnYKD9QSo0HnF65513sGXLFqxcuVLlQ3y9evXCmjVrMG3aNLXLCgoKwvjx42vN4+zsrFH7zp49C29vbwQGBuLNN9+sM7+vr6+0kVaLFi3QqlUr2Nra1rr7uVwu1+hLykRERE2Ca5x0TuPAaevWrfjkk08wePBgzJgxQ0p3dXXFH3/8oVFZVlZWsLKy0rQJNUpJSYG3tzcCAgIQHh6ucVsA4ODBg8jNzcWIESN01i4iIqImwcBJ5zQeg7t8+bL01eP7VVRUoKysTCeNqk5mZiaSkpKQmZmJ8vJyJCUlISkpCTdv3gRwL2gaNGgQfHx8EBISgpycHOTk5Khsm3D58mV069YNx48fl9IiIyNx7NgxpKWl4YsvvsDYsWMxd+5cdO3atcH6QkRERA8njUecevTogV9++QVOTk4q6V9//bX0teOGsGTJEmzZskU6r6zr0KFD8PLywtdff42rV69i27Zt2LZtm5TPyckJGRkZAICysjKkpqaiuLhYup6amoqFCxfi+vXrcHZ2xuLFizF37twG6wcREVGj4YiTzmn8rbrvvvsOkyZNwsKFC7Fs2TIsXboUqamp2Lp1K77//vtaF2Prq8pvAhVkpfGTK0T34ydXiKooLCyE0q5943yrbs5IKORafKuupAxt/reb36q7j8ZTdcOHD8eOHTuwd+9eyGQyLFmyBOfOncN33333rwyaiIiI6N+jXvs4DRkyBEOGDNF1W4iIiEiXOFWnc/UKnIiIiOghwMBJ59QKnFq3bq327uDXr1/XqkFEREREzZVagdPatWulP1+7dg3vvPMOhgwZAnd3dwD3Pm/y448/6vzDvkRERKQFjjjpnFqBU0BAgPTnMWPGYNmyZQgKCpLS5syZg/Xr1+PAgQN8lZ+IiKi50PazKfzkShUaP5Eff/wRQ4cOrZI+ZMgQHDhwQCeNIiIiImqONA6c2rRpg+jo6Crpu3fvRps2bXTSKCIiItKByqk6bQ5SofFbdUuXLsW0adMQFxcnrXE6duwYYmJi8Omnn+q8gURERFRPMmi5xklnLdEbGgdOU6ZMQffu3fG///0Pu3btghACLi4u+PXXX9G/f/+GaCMRERHVBxeH61y99nHq37+/yvfgiIiIiP4N1FrjVFhYqFGhRUVF9WoMERER6Y7MwEDrg1Sp9URat26N3NxctQtt164d/vzzz3o3ioiIiHRB24XhnKp7kFpTdUIIfPrppzA3N1er0LKyMq0aRURERNQcqRU4tW/fHhs3blS7UFtbW7Rs2bLejSIiIiId4OJwnVMrcMrIyGjgZhAREZHOMXDSOa76IiIiIlJTvbYjICIioocAv1Wnc3wiRERE+qqRP7myYcMGuLq6QqFQQKFQwN3dHfv27ZOuCyEQFhYGe3t7mJiYwMvLCykpKSpleHl5QSaTqRzjx4+vs+4PP/wQHTp0gLGxMfr164dffvlF5bo6dauDgRMRERHphIODA1asWIHExEQkJibC29sb/v7+UoCycuVKrF69GuvXr0dCQgJsbW3h4+NTZf/HwMBAZGdnS8fHH39ca707duxAcHAwFi9ejFOnTuHJJ5+En58fMjMzpTzq1l0XBk5ERET6qpFHnIYPH45hw4ahS5cu6NKlC8LDw2Fubo5jx45BCIG1a9di8eLFGD16NHr27IktW7aguLgY27dvVynH1NQUtra20qFUKmutd/Xq1Zg2bRpeeukldO/eHWvXroWjoyM2bNgAABrVXZd6BU6//PILXnjhBbi7u+Py5csAgM8//xyHDx+uT3FERETUEHQUOBUWFqocJSUldVZdXl6OqKgo3Lp1C+7u7khPT0dOTg58fX2lPHK5HJ6enjhy5IjKvdu2bYOVlRV69OiBefPm1ToqVFpaihMnTqiUCwC+vr5SuZrUXReNA6edO3diyJAhMDExwalTp6SHV1RUhIiICE2LIyIiooZSuThcmwOAo6MjlEqldCxfvrzGKpOTk2Fubg65XI4ZM2YgOjoaLi4uyMnJAQDY2Nio5LexsZGuAcDzzz+PL7/8EnFxcXjrrbewc+dOjB49usb68vLyUF5eXmu56tatDo3fqnvnnXfw0UcfYfLkyYiKipLSPTw8sGzZMk2LIyIiomYuKysLCoVCOpfL5TXm7dq1K5KSkpCfn4+dO3ciICAA8fHx0nXZA9N/QgiVtMDAQOnPPXv2ROfOneHm5oaTJ0+ib9++NdZbV7nq5qmLxiNOqampeOqpp6qkKxQK5Ofna1ocERERNRQdTdVVviVXedQWOBkZGaFTp05wc3PD8uXL0bt3b6xbtw62trYAUGWEJzc3t8pI0P369u2Lli1b4sKFC9Vet7KygqGhYa3l1rfu6mgcONnZ2eHixYtV0g8fPoxHHnlE0+KIiIiooTTy4vDqCCFQUlKCDh06wNbWFrGxsdK10tJSxMfHw8PDo8b7U1JSUFZWBjs7u2qvGxkZoV+/firlAkBsbKxUbn3rro7GU3XTp0/Hq6++is8++wwymQx///03jh49innz5mHJkiWaFkdERER6YtGiRfDz84OjoyOKiooQFRWFuLg4xMTEQCaTITg4GBEREejcuTM6d+6MiIgImJqaYuLEiQCAtLQ0bNu2DcOGDYOVlRXOnj2L1157DX369METTzwh1TN48GCMGjUKQUFBAICQkBBMmjQJbm5ucHd3xyeffILMzEzMmDEDANSqW10aB07z589HQUEBBg0ahDt37uCpp56CXC7HvHnzpA4QERFRM9DIO4dfuXIFkyZNQnZ2NpRKJVxdXRETEwMfHx8A92KI27dvY9asWbhx4wb69++P/fv3w8LCAsC90aOffvoJ69atw82bN+Ho6IhnnnkGoaGhMDQ0lOpJS0tDXl6edD5u3Dhcu3YNy5YtQ3Z2Nnr27Im9e/fCyclJylNX3eqSCSGERnf8o7i4GGfPnkVFRQVcXFxgbm5en2L0QmFhIZRKJQqy0qBQaPYXQKTXWhg1dQuImp3CwkIo7dqjoKBAZcG1zutQKnE9YjoUxjWvR6qznDslsFz0cYO29WFT72/VmZqaws3NTZdtISIiImrW1Aqcats/4UG7du2qd2OIiIhIh7Rd4K2DxeH6Rq3A6f6tzoUQiI6OhlKplEacTpw4gfz8fI0CLCIiImpgDJx0Tq3AKTIyUvrzggUL8Nxzz+Gjjz6SFmqVl5dj1qxZnP8kIiIivabxUvvPPvsM8+bNU1ndbmhoiJCQEHz22Wc6bRwRERFpQabl51ZkWryRp6c0fiJ3797FuXPnqqSfO3cOFRUVOmkUERER6UAz2ABT32j8Vt2LL76IqVOn4uLFixgwYAAA4NixY1ixYgVefPFFnTeQiIiI6olrnHRO48Dpvffeg62tLdasWYPs7GwA9z7DMn/+fLz22ms6byARERFRc6Fx4GRgYID58+dj/vz5KCwsBAAuCiciImqOZFquU+IapyrqvQEmwICJiIioWTOQ3Tu0uZ9UaBw4dejQAbJa5jz//PNPrRpERERE1FxpHDgFBwernJeVleHUqVOIiYnB66+/rqt2ERERkbY4VadzGgdOr776arXpH3zwARITE7VuEBEREekI36rTOZ2Fkn5+fti5c6euiiMiIiJqdrRaHH6/b775BpaWlroqjoiIiLRVuQO4NveTCo0Dpz59+qgsDhdCICcnB1evXsWHH36o08YRERGRFjhVp3MaB07+/v4qgZOBgQHatm0LLy8vdOvWTaeNIyIiImpONA6cwsLCGqAZREREpHN8q07nNH4ihoaGyM3NrZJ+7do1GBoa6qRRREREpAMyaPmR36buQPOj8YiTEKLa9JKSEhgZGWndICIiItIRLg7XObUDp//9738AAJlMhk8//RTm5ubStfLycvz8889c40RERER6Te3Aac2aNQDujTh99NFHKtNyRkZGcHZ2xkcffaT7FhIREVH98K06nVM7cEpPTwcADBo0CLt27ULr1q0brFFERESkA1wcrnMar3E6dOhQQ7SDiIiIqNlTK3AKCQnB22+/DTMzM4SEhNSad/Xq1TppGBEREWlJJgMMOFWnS2oFTqdOnUJZWRkA4OTJkyobYBIREVEzxak6nVMrcLp/ei4uLq6h2kJERETUrGkcSk6dOhVFRUVV0m/duoWpU6fqpFFERESkA1ptfqnlG3l6SuPAacuWLbh9+3aV9Nu3b2Pr1q06aRQRERHpQOVUnTYHqVD7rbrCwkIIISCEQFFREYyNjaVr5eXl2Lt3L6ytrRukkURERETNgdqhZKtWrWBpaQmZTIYuXbqgdevW0mFlZYWpU6di9uzZDdlWIiIi0oSBTPtDAxs2bICrqysUCgUUCgXc3d2xb98+6boQAmFhYbC3t4eJiQm8vLyQkpJSbVlCCPj5+UEmk2H37t211uvs7AyZTFbluD8umTJlSpXrAwYM0Kh/gAYjTocOHYIQAt7e3ti5cycsLS2la0ZGRnBycoK9vb3GDSAiIqIG0sg7hzs4OGDFihXo1KkTgHvLe/z9/XHq1Cn06NEDK1euxOrVq7F582Z06dIF77zzDnx8fJCamgoLCwuVstauXav2W/wJCQkoLy+Xzs+cOQMfHx+MHTtWJd/QoUMRGRkpndfnG7tqB06enp4A7u0g7ujoCAN++I+IiKh5a+TtCIYPH65yHh4ejg0bNuDYsWNwcXHB2rVrsXjxYowePRrAvcDKxsYG27dvx/Tp06X7Tp8+jdWrVyMhIQF2dnZ11tu2bVuV8xUrVqBjx45S7FJJLpfD1tZWoz49SOOn6eTkBAMDAxQXF+OPP/7A77//rnIQERGRfiksLFQ5SkpK6rynvLwcUVFRuHXrFtzd3ZGeno6cnBz4+vpKeeRyOTw9PXHkyBEprbi4GBMmTMD69evrFeSUlpbiiy++wNSpU6uMWMXFxcHa2hpdunRBYGAgcnNzNS5f48Dp6tWr+M9//gMLCwv06NEDffr0UTkaSnh4ODw8PGBqaopWrVpVuX7t2jUMHToU9vb2kMvlcHR0RFBQEAoLC2stt6SkBK+88gqsrKxgZmaGESNG4K+//mqgXhARETUiHa1xcnR0hFKplI7ly5fXWGVycjLMzc0hl8sxY8YMREdHw8XFBTk5OQAAGxsblfw2NjbSNQCYO3cuPDw84O/vX68u7969G/n5+ZgyZYpKup+fH7Zt24aDBw9i1apVSEhIgLe3t1pB4P00/lZdcHAwbty4gWPHjmHQoEGIjo7GlStX8M4772DVqlWaFqe20tJSjB07Fu7u7ti0aVOV6wYGBvD398c777yDtm3b4uLFi5g9ezauX7+O7du319qf7777DlFRUWjTpg1ee+01/Oc//8GJEydgaGjYYP0hIiJqcDKZllN19wKnrKwsKBQKKVkul9d4S9euXZGUlIT8/Hzs3LkTAQEBiI+Pv69I1VEgIYSUtmfPHhw8eBCnTp2qd5M3bdoEPz+/Kuuux40bJ/25Z8+ecHNzg5OTE3744Qdp6lAdGgdOBw8exLfffovHHnsMBgYGcHJygo+PDxQKBZYvX45nnnlG0yLVsnTpUgDA5s2bq73eunVrzJw5Uzp3cnLCrFmz8N///rfGMgsKCrBp0yZ8/vnnePrppwEAX3zxBRwdHXHgwAEMGTJEdx0gIiJ6SFW+JacOIyMjaXG4m5sbEhISsG7dOixYsAAAkJOTo7JuKTc3VxqFOnjwINLS0qrMLI0ZMwZPPvlknV8vuXTpEg4cOIBdu3bV2U47Ozs4OTnhwoULavWrksZh6K1bt6T9miwtLXH16lUAQK9evXDy5ElNi2swf//9N3bt2lVlYdj9Tpw4gbKyMpX5Vnt7e/Ts2VNlvvVBJSUlVeZ7iYiImp1msHO4EAIlJSXo0KEDbG1tERsbK10rLS1FfHw8PDw8AABvvPEGfv/9dyQlJUkHAKxZs0blbbiaREZGwtraWq1BnGvXriErK0utxef30zhw6tq1K1JTUwEAjz76KD7++GNcvnwZH330kcaVN4QJEybA1NQU7dq1g0KhwKefflpj3pycHBgZGaF169Yq6Q/Otz5o+fLlKnO9jo6OOms/ERGRzjTyzuGLFi3CL7/8goyMDCQnJ2Px4sWIi4vD888/D5lMhuDgYERERCA6OhpnzpzBlClTYGpqiokTJwIAbG1t0bNnT5UDANq3b48OHTpI9QwePBjr169XqbuiogKRkZEICAhAixaqE2o3b97EvHnzcPToUWRkZCAuLg7Dhw+HlZUVRo0apVEfNQ6cgoODkZ2dDQAIDQ1FTEwM2rdvj//973+IiIjQqKywsLBqN6y6/0hMTNSozDVr1uDkyZPYvXs30tLSEBISotH9gOp8a3UWLlyIgoIC6cjKytK4DiIiIn1z5coVTJo0CV27dsXgwYPx22+/ISYmBj4+PgCA+fPnIzg4GLNmzYKbmxsuX76M/fv3V9nDqS5paWnIy8tTSTtw4AAyMzOr/W6uoaEhkpOT4e/vjy5duiAgIABdunTB0aNHNa5bJoQQGt3xgMptCdq3bw8rKyuN7s3Ly6vS8Qc5OzurfN5l8+bNCA4ORn5+fp3lHz58GE8++ST+/vvvakfDDh48iMGDB+P69esqo069e/fGyJEjpXVVdSksLIRSqURBVhoUCs3+Aoj0WgvNN5cj0neFhYVQ2rVHQUGB2uuG6lWHUonrUaugMDWpfznFt2E5/rUGbevDRuPF4Q8yNTVF375963WvlZWVxsGWJipjwppeNezXrx9atmyJ2NhYPPfccwCA7OxsnDlzBitXrmywdhERETWKRt4A899ArcBJk+mu1atX17sxtcnMzMT169eRmZmJ8vJyacFYp06dYG5ujr179+LKlSt47LHHYG5ujrNnz2L+/Pl44okn4OzsDAC4fPkyBg8ejK1bt+Lxxx+HUqnEtGnT8Nprr6FNmzawtLTEvHnz0KtXL+ktOyIiIqJKagVO6u6noO43ZepjyZIl2LJli3ReudnmoUOH4OXlBRMTE2zcuBFz585FSUkJHB0dMXr0aLzxxhvSPWVlZUhNTUVxcbGUtmbNGrRo0QLPPfccbt++jcGDB2Pz5s3cw4mIiB5+jfytun8Drdc4Edc4EdWIa5yIqmjUNU7f/E/7NU7PzuEap/tovcaJiIiImitt92LiiNODuOqLiIiISE0ccSIiItJXfKtO5xg4ERER6SsuDtc5hpJEREREauKIExERkb4yMLh3aHM/qWDgREREpK84VadzDCWJiIiI1MQRJyIiIn0lk2n5Vh1HnB7EwImIiEhfcapO5zhVR0RERKQmjjgRERHpK26AqXMMnIiIiPSVgezeoc39pIKBExERkb7iiJPO8YkQERERqYkjTkRERPqKb9XpHAMnIiIifcWpOp3jEyEiIiJSE0eciIiI9JRMJoNMi+k2be7VVwyciIiI9BWn6nSOT4SIiIhITRxxIiIi0lcccdI5Bk5ERET6SqblzuFc41QFQ0kiIiIiNXHEiYiISF9xqk7nGDgRERHpK+4crnMMJYmIiPSVTPb/o071OjQLnDZs2ABXV1coFAooFAq4u7tj37590nUhBMLCwmBvbw8TExN4eXkhJSWl2rKEEPDz84NMJsPu3btrrTcsLEzas6rysLW1rVKeunXXhoETERER6YSDgwNWrFiBxMREJCYmwtvbG/7+/lKAsnLlSqxevRrr169HQkICbG1t4ePjg6KioiplrV27VqMNOHv06IHs7GzpSE5OVrmuSd21YeBERESkryqn6rQ5NDB8+HAMGzYMXbp0QZcuXRAeHg5zc3McO3YMQgisXbsWixcvxujRo9GzZ09s2bIFxcXF2L59u0o5p0+fxurVq/HZZ5+pXXeLFi1ga2srHW3btpWuaVJ3XRg4ERER6Sutpun+f2F5YWGhylFSUlJn1eXl5YiKisKtW7fg7u6O9PR05OTkwNfXV8ojl8vh6emJI0eOSGnFxcWYMGEC1q9fX2W6rTYXLlyAvb09OnTogPHjx+PPP/+UrqlbtzoYOBEREVGtHB0doVQqpWP58uU15k1OToa5uTnkcjlmzJiB6OhouLi4ICcnBwBgY2Ojkt/Gxka6BgBz586Fh4cH/P391W5f//79sXXrVvz444/YuHEjcnJy4OHhgWvXrgGA2nWrg2/VERER6SsDLTfA/OferKwsKBQKKVkul9d4S9euXZGUlIT8/Hzs3LkTAQEBiI+Pl64/uG5JCCGl7dmzBwcPHsSpU6c0aqafn5/05169esHd3R0dO3bEli1bEBISolbd6uKIExERkb7S0VRd5VtylUdtgZORkRE6deoENzc3LF++HL1798a6deukabcHR3hyc3OlkaCDBw8iLS0NrVq1QosWLdCixb3xnTFjxsDLy0vtbpuZmaFXr164cOECAKhVt7oYOBEREVGDEUKgpKQEHTp0gK2tLWJjY6VrpaWliI+Ph4eHBwDgjTfewO+//46kpCTpAIA1a9YgMjJS7TpLSkpw7tw52NnZAYBadauLU3VERET6qpE3wFy0aBH8/Pzg6OiIoqIiREVFIS4uDjExMZDJZAgODkZERAQ6d+6Mzp07IyIiAqamppg4cSIASG/EPah9+/bo0KGDdD548GCMGjUKQUFBAIB58+Zh+PDhaN++PXJzc/HOO++gsLAQAQEB/3Sj7rrVxcCJiIhIXzXyJ1euXLmCSZMmITs7G0qlEq6uroiJiYGPjw8AYP78+bh9+zZmzZqFGzduoH///ti/fz8sLCw0qictLQ15eXnS+V9//YUJEyYgLy8Pbdu2xYABA3Ds2DE4OTlJeXRVt0wIITS6g6ooLCyEUqlEQVYaFArN/gKI9FoLo6ZuAVGzU1hYCKVdexQUFKgsuNZ5HUolbhzdC4W5Wf3LuXkLrd2HNWhbHzYccSIiItJX/FadzjFwIiIi0leNPFX3b8DAiYiISF8ZGNw7tLmfVPCJEBEREamJI05ERER6SiaTabwz9oP3kyoGTkRERPpKJtNyjRMDpwdxqo6IiIhITRxxIiIi0lfcjkDnGDgRERHpLS23I+DEVBV8IkRERERq4ogTERGRvuJUnc4xcCIiItJX3ABT5xg4ERER6SuOOOkcQ0kiIiIiNXHEiYiISF/xI786x8CJiIhIX3GqTucYShIRERGpiSNOREREekv2z6HN/XQ/Bk5ERET6ilN1OsepOiIiIiI1ccSJiIhIX3HESecYOBEREektrnHSNU7VEREREamJI05ERET6ilN1OsfAiYiISF9xpk7nHpqpuvDwcHh4eMDU1BStWrWqcv3atWsYOnQo7O3tIZfL4ejoiKCgIBQWFtZa7ieffAIvLy8oFArIZDLk5+c3TAeIiIganUwHB93voQmcSktLMXbsWMycObPa6wYGBvD398eePXtw/vx5bN68GQcOHMCMGTNqLbe4uBhDhw7FokWLGqLZREREpEcemqm6pUuXAgA2b95c7fXWrVurBFVOTk6YNWsW/vvf/9ZabnBwMAAgLi5OF80kIiJqPrjGSecemsBJU3///Td27doFT0/Ppm4KERFR05BBy8BJZy3RGw/NVJ26JkyYAFNTU7Rr1w4KhQKffvqpzusoKSlBYWGhykFERET6r0kDp7CwMMhkslqPxMREjcpcs2YNTp48id27dyMtLQ0hISE6b/fy5cuhVCqlw9HRUed1EBERaa9xF4dv2LABrq6uUCgUUCgUcHd3x759+6TrQgiEhYXB3t4eJiYm8PLyQkpKSrVlCSHg5+cHmUyG3bt311rv8uXL8dhjj8HCwgLW1tYYOXIkUlNTVfJMmTKlSowxYMAAjfoHNPFUXVBQEMaPH19rHmdnZ43KtLW1ha2tLbp164Y2bdrgySefxFtvvQU7OzstWqpq4cKFKgFZYWEhgyciImp+GnmNk4ODA1asWIFOnToBALZs2QJ/f3+cOnUKPXr0wMqVK7F69Wps3rwZXbp0wTvvvAMfHx+kpqbCwsJCpay1a9dCpmb98fHxmD17Nh577DHcvXsXixcvhq+vL86ePQszMzMp39ChQxEZGSmdGxkZadQ/oIkDJysrK1hZWTVY+UIIAPem1nRJLpdDLpfrtEwiIqKH3fDhw1XOw8PDsWHDBhw7dgwuLi5Yu3YtFi9ejNGjRwO4F1jZ2Nhg+/btmD59unTf6dOnsXr1aiQkJKg18BETE6NyHhkZCWtra5w4cQJPPfWUlC6Xy2Fra6tNFx+eNU6ZmZlISkpCZmYmysvLkZSUhKSkJNy8eRMAsHfvXkRGRuLMmTPIyMjA3r17MXPmTDzxxBPSqNXly5fRrVs3HD9+XCo3JycHSUlJuHjxIgAgOTkZSUlJuH79eqP3kYiISLeabh+n8vJyREVF4datW3B3d0d6ejpycnLg6+sr5ZHL5fD09MSRI0ektOLiYkyYMAHr16+vd5BTUFAAALC0tFRJj4uLg7W1Nbp06YLAwEDk5uZqXPZD81bdkiVLsGXLFum8T58+AIBDhw7By8sLJiYm2LhxI+bOnYuSkhI4Ojpi9OjReOONN6R7ysrKkJqaiuLiYinto48+krY6ACBFppGRkZgyZUoD94qIiKgB6Wiq7sGXoGqbeUlOToa7uzvu3LkDc3NzREdHw8XFRQqObGxsVPLb2Njg0qVL0vncuXPh4eEBf3//ejVZCIGQkBAMHDgQPXv2lNL9/PwwduxYODk5IT09HW+99Ra8vb1x4sQJjWaRHprAafPmzTXu4QQAgwYNUolYq+Ps7CxN31UKCwtDWFiYDlpIRESknx5cxxsaGlrj786uXbsiKSkJ+fn52LlzJwICAhAfHy9df3DdkhBCStuzZw8OHjyIU6dO1butQUFB+P3333H48GGV9HHjxkl/7tmzJ9zc3ODk5IQffvhBmjpUx0MTOBEREZGGdDTilJWVBYVCISXXNkJjZGQkLQ53c3NDQkIC1q1bhwULFgC4t0Tm/nVLubm50ijUwYMHkZaWVuXTamPGjMGTTz5Z52bVr7zyCvbs2YOff/4ZDg4Otea1s7ODk5MTLly4UGu+Bz00a5yIiIhIU7pZ41S5vUDlocnUlhACJSUl6NChA2xtbREbGytdKy0tRXx8PDw8PAAAb7zxBn7//XdpHXNSUhKAe1sN3f82XHV1BAUFYdeuXTh48CA6dOhQZ7uuXbuGrKwsjd+654gTERGRnqrcr0ib+zWxaNEi+Pn5wdHREUVFRYiKikJcXBxiYmIgk8kQHByMiIgIdO7cGZ07d0ZERARMTU0xceJEAP+/pdCD2rdvrxIMDR48GKNGjUJQUBAAYPbs2di+fTu+/fZbWFhYICcnBwCgVCphYmKCmzdvIiwsDGPGjIGdnR0yMjKwaNEiWFlZYdSoURr1kYETERER6cSVK1cwadIkZGdnQ6lUwtXVFTExMfDx8QEAzJ8/H7dv38asWbNw48YN9O/fH/v376+yh1Nd0tLSkJeXJ51v2LABAODl5aWSr/JFL0NDQyQnJ2Pr1q3Iz8+HnZ0dBg0ahB07dmhct0w8uFqaNFZYWAilUomCrDQoFJr9BRDptRaaby5HpO8KCwuhtGuPgoIClXVDOq9DqURB+lkoNAwMVMopKoKyg0uDtvVhwxEnIiIivaXdXkz8ym9VXBxOREREpCaOOBEREektLbcj4IhTFQyciIiI9FUjf+T334BTdURERERq4ogTERGR3uLicF1j4ERERKSvOFWnc5yqIyIiIlITR5yIiIj0FWfqdI6BExERkd5i5KRrDJyIiIj0Fdc46RzXOBERERGpiSNORERE+oojTjrHwImIiEhvcY2TrnGqjoiIiEhNHHEiIiLSVzJoOVWns5boDQZORERE+oprnHSOU3VEREREauKIExERkd7i4nBdY+BERESkrzhVp3MMnHRACAEAKCwqauKWEDUzLVo2dQuImp3K3xWVvzsao66mul8fMXDSgaJ/frAcXR5t2oYQEdFDo6ioCEqlskHKNjIygq2tLRy79NC6LFtbWxgZGemgVfpBJhoj5NVzFRUV+Pvvv2FhYQGZng5rFhYWwtHREVlZWVAoFE3dnIcKn1398dnVH5+ddhry+QkhUFRUBHt7exgYNNw7Wnfu3EFpaanW5RgZGcHY2FgHLdIPHHHSAQMDAzg4ODR1MxqFQqHgf4Tric+u/vjs6o/PTjsN9fwaaqTpfsbGxgx4GgC3IyAiIiJSEwMnIiIiIjUxcCK1yOVyhIaGQi6XN3VTHjp8dvXHZ1d/fHba4fOjmnBxOBEREZGaOOJEREREpCYGTkRERERqYuBEREREpCYGTkRERERqYuBEkhEjRqB9+/YwNjaGnZ0dJk2ahL///lslT0JCAgYPHoxWrVqhdevW8PX1RVJSUo1lZmRkQCaTVXt8/fXXDdyjxtMQz67S0aNH4e3tDTMzM7Rq1QpeXl64fft2A/Wk8TXUs/Py8qryMzd+/PgG7Enja8ifO+DeDtd+fn6QyWTYvXu37jvQhBrq2U2fPh0dO3aEiYkJ2rZtC39/f/zxxx8N2BNqbAycSDJo0CB89dVXSE1Nxc6dO5GWloZnn31Wul5UVIQhQ4agffv2+O2333D48GEoFAoMGTIEZWVl1Zbp6OiI7OxslWPp0qUwMzODn59fY3WtwTXEswPuBU1Dhw6Fr68vjh8/joSEBAQFBTXoZxoaW0M9OwAIDAxU+dn7+OOPG7o7jaohnx0ArF27Vm8/I9VQz65fv36IjIzEuXPn8OOPP0IIAV9fX5SXlzdGt6gxCKIafPvtt0Imk4nS0lIhhBAJCQkCgMjMzJTy/P777wKAuHjxotrlPvroo2Lq1Kk6b29zoqtn179/f/Hmm282eHubE109O09PT/Hqq682dHObFV3+m01KShIODg4iOztbABDR0dEN2fQm11D/vTt9+rTG91Dzpj//s5V06vr169i2bRs8PDzQsmVLAEDXrl1hZWWFTZs2obS0FLdv38amTZvQo0cPODk5qVXuiRMnkJSUhGnTpjVk85uUrp5dbm4ufvvtN1hbW8PDwwM2Njbw9PTE4cOHG7M7jUrXP3fbtm2DlZUVevTogXnz5qGoqKgxutEkdPnsiouLMWHCBKxfvx62traN1YUm01D/vbt16xYiIyPRoUMHODo6NmQXqDE1deRGzcv8+fOFqampACAGDBgg8vLyVK6fOXNGdOzYURgYGAgDAwPRrVs3cenSJbXLnzlzpujevbuum90s6PrZHT16VAAQlpaW4rPPPhMnT54UwcHBwsjISJw/f76hu9OoGuLn7pNPPhGxsbEiOTlZfPnll8LZ2Vk8/fTTDdmNJtEQz+7ll18W06ZNk86hpyNODfXfuw8++ECYmZkJAKJbt24cbdIzDJz0XGhoqABQ65GQkCDlv3r1qkhNTRX79+8XTzzxhBg2bJioqKgQQghRXFwsHn/8cTF58mRx/PhxcfToUTFmzBjRo0cPUVxcXGdbiouLhVKpFO+9916D9VeXmvrZ/frrrwKAWLhwoUp6r169xBtvvNFwHdeBpn521UlMTBQAxIkTJ3TeX11q6mf37bffik6dOomioiIp7WEJnJr62VXKz88X58+fF/Hx8WL48OGib9++4vbt2w3ad2o8DJz03NWrV8W5c+dqPWr6B52VlSUAiCNHjgghhPj000+FtbW1KC8vl/KUlJQIU1NT8eWXX9bZlq1bt4qWLVuK3Nxc3XSugTX1s/vzzz8FAPH555+rpD/33HNi4sSJOuplw2jqZ1ediooK0bJlSxEVFaVd5xpYUz+7V199VchkMmFoaCgdAISBgYHw9PTUeX91qamfXXUq79m+fbt2naNmo4WupvyoebKysoKVlVW97hX/fMawpKQEwL11DwYGBipv2VSeV1RU1Fnepk2bMGLECLRt27Ze7WlsTf3snJ2dYW9vj9TUVJX08+fPN/s3Epv62VUnJSUFZWVlsLOzq1e7GktTP7s33ngDL730kkpar169sGbNGgwfPrxe7WosTf3saiu7slzSA00ZtVHz8dtvv4n3339fnDp1SmRkZIiDBw+KgQMHio4dO4o7d+4IIYQ4d+6ckMvlYubMmeLs2bPizJkz4oUXXhBKpVL8/fffQggh/vrrL9G1a1fx22+/qZR/4cIFIZPJxL59+xq9bw2tIZ/dmjVrhEKhEF9//bW4cOGCePPNN4WxsbHerJloqGd38eJFsXTpUpGQkCDS09PFDz/8ILp16yb69Okj7t6922T91aWG/jd7PzwkU3Xqaqhnl5aWJiIiIkRiYqK4dOmSOHLkiPD39xeWlpbiypUrTdZf0i0GTiSEuPea7aBBg4SlpaWQy+XC2dlZzJgxQ/z1118q+SrXAiiVStG6dWvh7e0tjh49Kl1PT08XAMShQ4dU7lu4cKFwcHBQGfbWFw397JYvXy4cHByEqampcHd3F7/88ktjdKtRNNSzy8zMFE899ZSwtLQURkZGomPHjmLOnDni2rVrjdm9BtXQP3f307fAqaGe3eXLl4Wfn5+wtrYWLVu2FA4ODmLixInijz/+aMzuUQOTCfHP+CQRERER1Yr7OBERERGpiYETERERkZoYOBERERGpiYETERERkZoYOBERERGpiYETERERkZoYOBERERGpiYET0T+8vLwQHBysV/VOmTIFI0eO1KoMZ2dnyGQyyGQy5Ofn15hv8+bNaNWqlVZ1Uc2mTJki/T3s3r27qZtD9K/FwImoie3atQtvv/22dO7s7Iy1a9c2XYOqsWzZMmRnZ0OpVDZ1U/ReXFxctUHqunXrkJ2d3TSNIiIJP/JL1MQsLS2bugl1srCwgK2tbVM3AwBQVlaGli1bNnUzGp1SqWTgStQMcMSJqAY3btzA5MmT0bp1a5iamsLPzw8XLlyQrldOTf3444/o3r07zM3NMXToUJVRgbt372LOnDlo1aoV2rRpgwULFiAgIEBl+uz+qTovLy9cunQJc+fOlaZlACAsLAyPPvqoSvvWrl0LZ2dn6by8vBwhISFSXfPnz8eDX1QSQmDlypV45JFHYGJigt69e+Obb76p1/PZvHkz2rdvD1NTU4waNQrXrl2rkue7775Dv379YGxsjEceeQRLly7F3bt3pet//PEHBg4cCGNjY7i4uODAgQMqU1EZGRmQyWT46quv4OXlBWNjY3zxxRcAgMjISHTv3h3Gxsbo1q0bPvzwQ5W6L1++jHHjxqF169Zo06YN/P39kZGRIV2Pi4vD448/DjMzM7Rq1QpPPPEELl26pFbf6+rX6tWr0atXL5iZmcHR0RGzZs3CzZs3peuXLl3C8OHD0bp1a5iZmaFHjx7Yu3cvMjIyMGjQIABA69atIZPJMGXKFLXaRESNg4ETUQ2mTJmCxMRE7NmzB0ePHoUQAsOGDUNZWZmUp7i4GO+99x4+//xz/Pzzz8jMzMS8efOk6++++y62bduGyMhI/PrrrygsLKx1fcquXbvg4OAgTY1pMjWzatUqfPbZZ9i0aRMOHz6M69evIzo6WiXPm2++icjISGzYsAEpKSmYO3cuXnjhBcTHx6v/YAD89ttvmDp1KmbNmoWkpCQMGjQI77zzjkqeH3/8ES+88ALmzJmDs2fP4uOPP8bmzZsRHh4OAKioqMDIkSNhamqK3377DZ988gkWL15cbX0LFizAnDlzcO7cOQwZMgQbN27E4sWLER4ejnPnziEiIgJvvfUWtmzZAuDe38ugQYNgbm6On3/+GYcPH5YC29LSUty9excjR46Ep6cnfv/9dxw9ehQvv/yyFKjWpq5+AYCBgQH+97//4cyZM9iyZQsOHjyI+fPnS9dnz56NkpIS/Pzzz0hOTsa7774Lc3NzODo6YufOnQCA1NRUZGdnY926dRr93RBRA2vSTwwTNSOenp7i1VdfFUIIcf78eQFA/Prrr9L1vLw8YWJiIr766ishhBCRkZECgLh48aKU54MPPhA2NjbSuY2Njfjvf/8rnd+9e1e0b99e+Pv7V1uvEEI4OTmJNWvWqLQtNDRU9O7dWyVtzZo1wsnJSTq3s7MTK1askM7LysqEg4ODVNfNmzeFsbGxOHLkiEo506ZNExMmTKjxuVTXngkTJoihQ4eqpI0bN04olUrp/MknnxQREREqeT7//HNhZ2cnhBBi3759okWLFiI7O1u6HhsbKwCI6OhoIcT/f31+7dq1KuU4OjqK7du3q6S9/fbbwt3dXQghxKZNm0TXrl1FRUWFdL2kpESYmJiIH3/8UVy7dk0AEHFxcTX2uyZ19as6X331lWjTpo103qtXLxEWFlZt3kOHDgkA4saNG9Vev//5EFHj4xonomqcO3cOLVq0QP/+/aW0Nm3aoGvXrjh37pyUZmpqio4dO0rndnZ2yM3NBQAUFBTgypUrePzxx6XrhoaG6NevHyoqKnTa3oKCAmRnZ8Pd3V1Ka9GiBdzc3KTpurNnz+LOnTvw8fFRube0tBR9+vTRqL5z585h1KhRKmnu7u6IiYmRzk+cOIGEhASVkZjy8nLcuXMHxcXFSE1NhaOjo8raqfuf1f3c3NykP1+9ehVZWVmYNm0aAgMDpfS7d+9Ka4BOnDiBixcvwsLCQqWcO3fuIC0tDb6+vpgyZQqGDBkCHx8fPP3003juuedgZ2dXZ9/r6pepqSkOHTqEiIgInD17FoWFhbh79y7u3LmDW7duwczMDHPmzMHMmTOxf/9+PP300xgzZgxcXV3rrJuImh4DJ6JqiAfWBt2ffv90zoOLlGUyWZV7H5z+qans2hgYGFS57/4pQ3VUBms//PAD2rVrp3JNLpdrVJY6faioqMDSpUsxevToKteMjY2rPMvamJmZqZQLABs3blQJbIF7gWllnn79+mHbtm1Vymrbti2Ae2uk5syZg5iYGOzYsQNvvvkmYmNjMWDAAK36denSJQwbNgwzZszA22+/DUtLSxw+fBjTpk2T/s5eeuklDBkyBD/88AP279+P5cuXY9WqVXjllVfUeh5E1HQYOBFVw8XFBXfv3sVvv/0GDw8PAMC1a9dw/vx5dO/eXa0ylEolbGxscPz4cTz55JMA7o1MnDp1qspC7/sZGRmhvLxcJa1t27bIyclRCTaSkpJU6rKzs8OxY8fw1FNPAbg3AnPixAn07dtX6pNcLkdmZiY8PT3V6kNNXFxccOzYMZW0B8/79u2L1NRUdOrUqdoyunXrhszMTFy5cgU2NjYAgISEhDrrtrGxQbt27fDnn3/i+eefrzZP3759sWPHDlhbW0OhUNRYVp8+fdCnTx8sXLgQ7u7u2L59e52BU139SkxMxN27d7Fq1SoYGNxbRvrVV19Vyefo6IgZM2ZgxowZWLhwITZu3IhXXnkFRkZGAFDlZ4CImgcGTkTV6Ny5M/z9/REYGIiPP/4YFhYWeOONN9CuXTv4+/urXc4rr7yC5cuXo1OnTujWrRvef/993Lhxo9aRFmdnZ/z8888YP3485HI5rKys4OXlhatXr2LlypV49tlnERMTg3379qkEBa+++ipWrFiBzp07o3v37li9erXKXkAWFhaYN28e5s6di4qKCgwcOBCFhYU4cuQIzM3NERAQoHa/5syZAw8PD6xcuRIjR47E/v37VabpAGDJkiX4z3/+A0dHR4wdOxYGBgb4/fffkZycjHfeeQc+Pj7o2LEjAgICsHLlShQVFUmLw+saiQoLC8OcOXOgUCjg5+eHkpISJCYm4saNGwgJCcHzzz+P//73v/D398eyZcvg4OCAzMxM7Nq1C6+//jrKysrwySefYMSIEbC3t0dqairOnz+PyZMn19n3uvrVsWNH3L17F++//z6GDx+OX3/9FR999JFKGcHBwfDz80OXLl1w48YNHDx4UArInZycIJPJ8P3332PYsGEwMTGBubm52n83RNTAmmx1FVEz8+Ai7evXr4tJkyYJpVIpTExMxJAhQ8T58+el65GRkSqLoYUQIjo6Wtz/z6qsrEwEBQUJhUIhWrduLRYsWCDGjh0rxo8fX2O9R48eFa6urkIul6uUtWHDBuHo6CjMzMzE5MmTRXh4uMri8LKyMvHqq68KhUIhWrVqJUJCQsTkyZNVFqJXVFSIdevWia5du4qWLVuKtm3biiFDhoj4+Pgan0t1i8OFuLcA28HBQZiYmIjhw4eL9957r8rziImJER4eHsLExEQoFArx+OOPi08++US6fu7cOfHEE08IIyMj0a1bN/Hdd98JACImJkYI8f+Lw0+dOlWl/m3btolHH31UGBkZidatW4unnnpK7Nq1S7qenZ0tJk+eLKysrIRcLhePPPKICAwMFAUFBSInJ0eMHDlS2NnZCSMjI+Hk5CSWLFkiysvLa3wOmvRr9erVws7OTvq52bp1q8qC76CgINGxY0chl8tF27ZtxaRJk0ReXp50/7Jly4Stra2QyWQiICBApW5wcThRk5IJUY8FF0RULxUVFejevTuee+45ld3CmzNnZ2cEBwc3yudofv31VwwcOBAXL15UWXRP/08mkyE6OlrrT+kQUf1wHyeiBnTp0iVs3LgR58+fR3JyMmbOnIn09HRMnDixqZumkQULFsDc3BwFBQU6LTc6OhqxsbHIyMjAgQMH8PLLL+OJJ55g0FSNGTNmcMqOqBngiBNRA8rKysL48eNx5swZCCHQs2dPrFixQlrA/TC4dOmS9DbYI488Ii141oWtW7fi7bffRlZWFqysrPD0009j1apVaNOmjc7q0FSPHj1q3EH8448/rnFBekPLzc1FYWEhgHvbXtz/piERNR4GTkRE97k/UHyQjY1Nlb2hiOjfhYETERERkZq4xomIiIhITQyciIiIiNTEwImIiIhITQyciIiIiNTEwImIiIhITQyciIiIiNTEwImIiIhITQyciIiIiNT0f/c2ATvrPb4nAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Find the time step with the maximum temperature\n", + "max_temp_index = era_5_hottest_day.sel(variable='t2m').argmax(dim='valid_time')\n", + "\n", + "# # Select the data corresponding to the maximum time and plot\n", + "era_5_hottest_day.sel(variable='t2m').isel(valid_time=max_temp_index).plot(cmap='Reds')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Save to file" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + " \n", + " " + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Create a data folder if it does not exist\n", + "if not os.path.exists('data'):\n", + " os.makedirs('data')" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "\n", + " \n", + " " + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Save the era_5_hottest_day temperature to a tif file\n", + "era_5_hottest_day.sel(variable='t2m').rio.to_raster(raster_path='data/era_5_hottest_day_t2m.tif', driver=\"COG\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "cities-cif", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.14" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/setup.py b/setup.py index 5fec8d6..a46fcde 100644 --- a/setup.py +++ b/setup.py @@ -23,7 +23,8 @@ "dask>=2023.11.0", "boto3", "overturemaps", - "scikit-learn>=1.5.1", + "cdsapi", + "timezonefinder", "scikit-image>=0.24.0", "exactextract>=0.2.0" ], diff --git a/tests/test_layers.py b/tests/test_layers.py index d26a917..f8e1fe2 100644 --- a/tests/test_layers.py +++ b/tests/test_layers.py @@ -5,6 +5,7 @@ AlosDSM, AverageNetBuildingHeight, BuiltUpHeight, + Era5HottestDay, EsaWorldCover, EsaWorldCoverClass, HighLandSurfaceTemperature, @@ -47,6 +48,11 @@ def test_built_up_height(): data = BuiltUpHeight().get_data(BBOX) assert np.size(data) > 0 +@pytest.mark.skip(reason="CDS API needs personal access token file to run") +def test_era_5_hottest_day(): + data = Era5HottestDay().get_data(BBOX) + assert np.size(data) > 0 + def test_esa_world_cover(): land_cover_class = EsaWorldCoverClass.BUILT_UP data = EsaWorldCover(land_cover_class=land_cover_class).get_data(BBOX) diff --git a/tests/test_metrics.py b/tests/test_metrics.py index dcd731b..8fd42cc 100644 --- a/tests/test_metrics.py +++ b/tests/test_metrics.py @@ -1,5 +1,6 @@ from city_metrix import * from .conftest import ZONES +import pytest def test_built_land_with_high_lst(): @@ -23,6 +24,12 @@ def test_built_land_without_tree_cover(): assert expected_zone_size == actual_indicator_size +@pytest.mark.skip(reason="CDS API needs personal access token file to run") +def test_era_5_met_preprocess(): + indicator = era_5_met_preprocessing(ZONES) + assert len(indicator) == 24 + + def test_mean_tree_cover(): indicator = mean_tree_cover(ZONES) expected_zone_size = ZONES.geometry.size