diff --git a/.DS_Store b/.DS_Store new file mode 100644 index 0000000..1aa90d7 Binary files /dev/null and b/.DS_Store differ diff --git a/data_access/.DS_Store b/data_access/.DS_Store new file mode 100644 index 0000000..21199e6 Binary files /dev/null and b/data_access/.DS_Store differ diff --git a/data_access/2005-06-08_00_00_00_2006-01-06_23_59_59-SST14NA-anal_temp-15_60_220_310-.5.NETCDF b/data_access/2005-06-08_00_00_00_2006-01-06_23_59_59-SST14NA-anal_temp-15_60_220_310-.5.NETCDF new file mode 100644 index 0000000..e7decdf Binary files /dev/null and b/data_access/2005-06-08_00_00_00_2006-01-06_23_59_59-SST14NA-anal_temp-15_60_220_310-.5.NETCDF differ diff --git a/data_access/README_noaa_class_sst.md b/data_access/README_noaa_class_sst.md new file mode 100644 index 0000000..639330c --- /dev/null +++ b/data_access/README_noaa_class_sst.md @@ -0,0 +1,63 @@ +**Author:** Johanna Kasischke +# LAB BOOK for homework 1 +**Course:** ESDP +**Due Date:** 04.12.2025 + +## Homework assignment 1: Overview +The objective for this homework is to write a python code for downloading large volumes of data from environmental data servers. In more detail, the task requires to develop a storage concept and strategy for managing the downloads, including keeping track of which files have been transferred, as well as identifying those that failed. + +## Selection of data server +The chosen variable for this exercise is the sea surface temperature (SST), which has a critical impact on hurricane development. Since the earth is warming, the average sea surface temperature also rises, which may lead to more or heavier development of hurricanes. The focus area is the atlantic. Different servers were reviewed to find those that provided the necessary data. After comparing several options, NOAA CLASS satellite data was selected. The final dataset used comes from the NESDIS satellite programme, which provides historical SST data. + +### Initial attempt: NOAA CLASS GOES SST data +The intended datasat to use for this exercise was the Coastal Geostationary SST products from NOAA GOES and Japanese MTSAT-1R satellites (https://www.class.noaa.gov/saa/products/search?sub_id=0&datatype_family=GOESSST&submit.x=19&submit.y=4). The advantage of this dataset was the high temporal resolution. The focus time period should have been the recent hurricane Melissa in the carribean. Since the extraction of the received data didn't work out, a different dataset from the NOAA CLASS website was used. Besides, the download process for the GOES data will be shortly explained for other user to try their best in extracting the data. +It is first necessary to specify a time period. The advanced search function is not useful, because whenever the user chooses one of the "Datatype" or "Region" options, no data is found. This is suboptimal, because data must be downloaded, which may not be necessary. After defining a time period, the search button can be selected. A new window will then appear in which the datasets that are relevant to the user can be selected. In order to download these datasets, it is first necessary to register and then log in to one's account. Subsequently, the data can be requested for download. Upon completion of the process, the user is sent an e-mail a few hours later, in which it is stated that the data can be downloaded via FTP or via a link. After downloading the recipient will receive an ".gz" file, which requires decompression. It is unfortunate, but the file contained within the folder is not of a specific file format, and thus cannot be read by Linux programmes or Python scripts. Because of that the used dataset was changed for demonstrating purposes of this exercise. + +### Final choice: NESDIS SST14NA dataset +Using this data, a closer look at how sea surface temperatures influence hurricane development in the Atlantic was taken. As the new data source only provided data up to 2016, the time period of interest was changed to the 2005 Atlantic Hurricane Season. This season lasted from 8 June 2005 to 6 January 2006. This season is historically significant (e.g. Hurricanes Katrina, Rita and Wilma). +A note on the website says that the product has been discontinued and that the alternative product is: Geo-Polar 5 km Global Night-time Only Blended SST. A big advantage of this dataset was, that the format netCDF was available and compatible with scientific python libraries (e.g. xarray). The used dataset can be also found in the folder for easy access. + +Please note: Adjusting the time period based on dataset availability is not ideal in real research. However, for this exercise, the workflow remains valid and transferable to other datasets. + +## Download of the specific data +The NESDIS SST14NA dataset was accessed via: https://www.aev.class.noaa.gov/saa/products/welcome. The following steps need to be followed to download the data: +1. Navigate to the portal and select "Sea surface temperature (14 km North America)" either at the top search bar or on the right-hand side of the website. +2. Specify the temporal range (here: June 8 2005 till January 6 2006). +3. Choose output format (NetCDF preferred). Other options include ASCII, GIF, spreadsheet, tab-delimited or CSV. +4. Select the distributed variable (SST). +5. Initiate the search --> the server retrieves the specific data +6. The download starts by itself. + +## Storage and management concept +To handle large datasets efficiently, the FAIR principles should be applied. Therefore it is necessary to design a useful directory hierarchy, in which raw and processed data are stored apart from each other. In a log-file changes to the data can be written down. +A metadata file should include all important information of the downloaded datafiles, such as download date, the status of the download, a source url and other necessary information. Additionally, an error handling strategy should be implemented. There should be a workflow about how to deal with failed downloads, corrupted files and duplicates. + +## Evaluation of the respective data portal +Finding the appropiate dataset was uncomplicated, but the original data file caused technical difficulties. As a result, the dataset has to be switched, which is not ideal, because the new dataset provdies fewer data points and lacks the detail required for a real analysis. While general technical details about the dataset were easily available, obtaining more comprehensive or in-depth documentation was not straightforward. Another weakness of the dataset might be the discontinuation, especially, if one would want to use data which is younger than 2016. + +## Scientific context +Why does the sea surface temperature matters for hurricanes? +For a hurricane to develop, the SST needs to be greater than 26 °C according to Smith, R. K. and Montgomery, M. T. (2023). Because of the high temperatures in the carribean, tropical storms are very likely to develop there. An increase in SST would probably intensify the development of hurricanes in this region. The hurricane season of 2005 was exceptionally active, which makes it a strong case study. + +There are some limitations of the SST14NA dataset. The regional coverage might exclude some parts of the atlantic, which might be necessary for a proper analysis. The resolution is only 14 km, so fine-scale ocean features might be missed. Because of the discontinued dataset, a comparison with more recent events is not possible. + +# How to scale data access and make code more user-friendly +For this homework, only a limited dataset (one season) was downloaded. In real-world applications, researchers may need multi-year or global datasets, which introduces challenges: +1. Scalability: + The python scripts (or any other programming language) must handle large volumes of data efficiently. Therefore, parallel downloads or an request prompt inside the python file could speed up the access. +2. Storage: + Large datasets require structured storage solutions such as large servers at the office or cloud storage. + Access more data is with this dataset not a problem. With the GOES dataset, more datafiles would have been needed to download, which would have caused greater file sizes. If the file sizes would have been greater, there are several things to do to organize the data. This includes: +- group the data by year, month, day etc. +- group the data by region +- group the data by specific events inside the time period +4. user-friendly code: + The provided scripts should incluede clear instructions in how to access the data and change the download volume. Therefore a good documentation is needed. The script itself should kept clear and simple. +5. Limits: + There may be some limits to downloading the data from a source. This can be download quotas or server restrictions, network bandwidth or local storage capacity. + +Applied to the used dataset, no information could be found in restrictions according to the maximum download quotas. + +# Literature +Smith, R. K. and Montgomery, M. T. (2023). Chapter 1 - observations of tropical cyclones. In Smith, R. K. and Montgomery, M. T., editors, Tropical Cyclones, volume 4 of Developments in Weather and Climate Science, pages 1–34. Elsevier. + diff --git a/data_access/homework2/.DS_Store b/data_access/homework2/.DS_Store new file mode 100644 index 0000000..60d3838 Binary files /dev/null and b/data_access/homework2/.DS_Store differ diff --git a/data_access/homework2/ERA5_control_flow.py b/data_access/homework2/ERA5_control_flow.py new file mode 100644 index 0000000..4efb742 --- /dev/null +++ b/data_access/homework2/ERA5_control_flow.py @@ -0,0 +1,124 @@ +import cdsapi +from datetime import datetime, timedelta +import healpy as hp +import numpy as np +import xarray as xr +import os +import zarr + + +def retrieve_and_process(dataset, start, end, variable, time, pressure_level, + data_format="netcdf", download_format = "unarchived"): + + # SETUP + START_DATE = datetime.strptime(start, "%Y-%m-%d") + END_DATE = datetime.strptime(end, "%Y-%m-%d") + c = cdsapi.Client() + current = START_DATE + + # DAILY LOOP + while current <= END_DATE: + date_str = current.strftime("%Y-%m-%d") + print(f"\nProcessing: {date_str}") + + # DOWNLOAD + filename = f"era5_{date_str}.nc" + c.retrieve(dataset, { + "product_type": ["reanalysis"], + "variable": variable, + "date": date_str, + "time": time, + "pressure_level": pressure_level, + "data_format": data_format, + "download_format": download_format + }).download(filename) + + + # INTERPOLATE DATA TO HEALPIX GRID + + ## Load datafile + ds = xr.open_dataset(filename) + + ## Get time and coordinates, rename to standard names + if 'valid_time' in ds.coords: + time_coord = ds['valid_time'].rename({'valid_time': 'time'}) + else: + time_coord = ds['time'] + + if 'pressure_level' in ds.coords: + level_coord = ds['pressure_level'].rename({'pressure_level': 'level'}) + elif 'level' in ds.coords: + level_coord = ds['level'] + else: + level_coord = None + + ## Implement and convert to healpix grid + lats, lons = ds['latitude'].values, ds['longitude'].values + lon_grid, lat_grid = np.meshgrid(lons, lats) + theta = np.radians(90.0 - lat_grid).flatten() + phi = np.radians(lon_grid).flatten() + + # Interpolate to both NSIDE values + for nside in [8, 16]: + print(f" NSIDE={nside}") + pix_indices = hp.ang2pix(nside, theta, phi) + npix = hp.nside2npix(nside) + zarr_path = f"era5_nside{nside}.zarr" # create path to store the zarr files + + ### loop through variables (here: q for specific humidity) + data_vars = {} + for var_name in ds.data_vars: + if var_name in ['latitude', 'longitude', 'number', 'expver']: + continue + + data = ds[var_name].values + healpix_map = np.full(list(data.shape[:-2]) + [npix], hp.UNSEEN, dtype=np.float32) + data_flat = data.reshape(-1, data.shape[-2] * data.shape[-1]) + healpix_flat = healpix_map.reshape(-1, npix) + + for i in range(data_flat.shape[0]): + valid = ~np.isnan(data_flat[i]) + if np.any(valid): + counts = np.bincount(pix_indices[valid], minlength=npix) + sums = np.bincount(pix_indices[valid], weights=data_flat[i, valid], minlength=npix) + healpix_flat[i, counts > 0] = sums[counts > 0] / counts[counts > 0] + + # Create DataArray to store data and metadata + if level_coord is not None and len(level_coord) > 1: + data_vars[var_name] = xr.DataArray( + healpix_flat.reshape(healpix_map.shape), + dims=['time', 'level', 'healpix'], + coords={'time': time_coord, 'level': level_coord, 'healpix': np.arange(npix)} + ) + else: + data_vars[var_name] = xr.DataArray( + healpix_flat.reshape(healpix_map.shape), + dims=['time', 'healpix'], + coords={'time': time_coord, 'healpix': np.arange(npix)} + ) + + # create xarray dataset + ds_hp = xr.Dataset(data_vars) + ds_hp.attrs['nside'] = nside + + # Save to zarr + if os.path.exists(zarr_path): + # For appending, use append_dim and reconsolidate after + ds_hp.to_zarr(zarr_path, mode='a', append_dim='time') + # Reconsolidate metadata after appending + zarr.consolidate_metadata(zarr_path) + else: + # For new files, set up chunking and consolidated metadata + encoding = {} + for var in data_vars: + if 'level' in ds_hp[var].dims: + encoding[var] = {'chunks': (1, len(level_coord), npix)} + else: + encoding[var] = {'chunks': (1, npix)} + + ds_hp.to_zarr(zarr_path, mode='w', encoding=encoding, consolidated=True) + + print(f" Saved to {zarr_path}") + + ds.close() + current += timedelta(days=1) \ No newline at end of file diff --git a/data_access/homework2/README_hw2_kasischke.md b/data_access/homework2/README_hw2_kasischke.md new file mode 100644 index 0000000..c96bdd8 --- /dev/null +++ b/data_access/homework2/README_hw2_kasischke.md @@ -0,0 +1,114 @@ +# Homework 2 +Author: Johanna Kasischke +Module: Earth System Data Processing + +## Part 1: Set-up of control flow +A control flow determines which part of the program, e.g. functions or statements, executes next. Some programming languages contain so called "control statements" (e.g. if/else statements). By default, a program runs line by line from top to bottom. With control structures (conditionals, loops, function calls), the default flow is altered. + +### Objective +For the first part of the homework, we need to set-up a control flow that handles daily data. +The workflow will be tested with a mock processing that doesn't download any data. + +### Code implementation +The mock workflow can be executed by running the file `control_flow.py`. + +At the beginning of the script the following modules are imported: `argparse`, `random`, and the `datetime, timedelta, date` subclasses from the `datetime` module. +The only parameter that must be defined initially is the start date of the project (e.g. `date(2026, 1, 1)`, which represents the beginning of the processing period. +Several mock functions are then defined to represent a realistic data-processing pipeline. +The functions are: `download_data`, `process_data`, `archive_data` and `run_pipeline` +These functions do not perform real data operations but are used to test the logic of the workflow. + +**The `main()` function** + +The function `main` executes the code above in the following way: + +1. The Argument Parser +With the `main()`function, the script can be runned from the command line with a specific date. The specific date can be included as a string and the lambda function converts the input string automatically into a Python date object. + +There are two sencarios that are handled: +1. Call script with date argument +If a date argument is provided, the script just executes `run_pipeline()` for the specific day. +2. No data provided +If no date is provided, the script is looking for missing data until today and runs the pipeline for every missing day. If the `run_pipeline()` function fails (randomly) for a day, the entire loop breaks. + +At the end of the script, two datasets were already added to the ARCHIVE folder to simulate a real scenario, where the code needs to look for missing data files. + +## Part 2: Download ERA5 humidity data +For the next part of the exercise, ERA5 humidity data is downloaded for the period 2024-12-01 to 2024-12-05 in 6-hourly intervals on pressure levels 975, 900, 800, 500, 300 hPa in the original lat-lon resolution. + +**Dataset description** + +For this task the **ERA5 hourly data on pressure levels from 1940 to present** dataset is used. + +The “ERA5 hourly data on pressure levels from 1940 to present" dataset is part of the fifth generation reanalysis provided by the ECMWF. The ERA5 dataset offers hourly estimates for a large number of atmospheric, ocean-wave and land-surface quantities from 1940 until today. This dataset is regridded to a regular lat-lon of 0.25 degrees for the reanalysis. It has a vertical coverage from 1000 hPa to 1hPa and the vertical resolution is 37 pressure levels. The temporal resolution is hourly. The data can be found in the Copernicus data store or under the following link: [https://cds.climate.copernicus.eu/datasets/reanalysis-era5-pressure-levels?tab=overview](https://cds.climate.copernicus.eu/datasets/reanalysis-era5-pressure-levels?tab=overview) + +### Code structure: + + +1. `load_ERA5_data.py` script +Here, the following parameters can be changed: + - variable (e.g. specific humidity) + - pressure level + - start and end date of the download + - if needed, the data_format and download_format + +All other parameters remain fixed in the `ERA5_control_flow.py` file. + +2. `ERA5_control_flow.py` script + +The script contains the main function `retrieve_and_process()`, which performs the complete workflow of downloading, processing, interpolation and saving the data. +To work with the script, the following modules need to be installed: + - `cdsapi` for downloading the ERA5 data + - `datetime`for handling dates + - `healpy`for interpolating the data into a healpix grid + - `numpy` for working with arrays + - `xarray` for working with netCDF files + - `os`for file operations + - `zarr` for using zarr as a storage method + +#### Daily Processing Workflow + +A daily loop is implemented to process one day at a time. For each day: + +1. ERA5 data is requested from the CDS API for the specific day, with chosen variables, times and pressure levels +2. The current netCDF file is loaded with xarray +3. The coordinates are converted into HEALPix coordinates: + - theta = colatitude (90° - latitude) in radians + - phi = longitude in radians + - `.flatten()´converts the 2D grids to 1D arrays +4. For each HEALPix resolution, an empty HEALPix array (time × level × npix) is created and filled with `hp.UNSEEN` + +#### Interpolation to HEALPix Grid +The interpolation is performed seperately for each time and pressure-level combination. +- First it is searched for valid (non-NaN) data points. + - `np.bincount` with pix_indices: counts how many grid points map to each HEALPix pixel + - `np.bincount` with weights: sums the values for each HEALPix pixel +- Then the sum is divided by count to get the average value per pixel. +The interpolated data is reshaped back into their original dimensions and stored in an xarray DataArray. + +#### Zarr Output and Chunking strategy +After interpolation, all variables are combined into a single xarray dataset and the HEALPix resolution (NSIDE) is stores as metadata. + +Then the data is saved in Zarr format using the following strategy: +- Chunking + - One day per chunk along the `time` dimension + - All pressure levels and HEALPix pixels are stored together +- Appending + - If file exists: append new day's data along time dimension + - If new file: create with chunking strategy (1 day per chunk) and consolidated metadata +- Metadata + - Consolidated metadata allows xarray to quickly read dimensions without loading all data + + +#### HEALPix Coordinate Convention +The standard coordinates are the colatitude (theta),0 at the North Pole, pi/2 at the equator and pi at the South Pole and the longitude (Phi) between 0 and 2 pi eastward, in a Mollview projection,phi=0 is at the center and increases eastward toward the left of the map. When visualized using a Mollview projection (`hp.mollview`), phi = 0 is located at the center of the map and increases eastward toward the left. + +## Comment on the plots +The resolution of the regridded data is not sufficient. I think, I might have an issue with the coordinates and the HealPix grid. But the computation time for creating the plots is much faster than with the original netCDF files. + + +## Tools +Claude AI was used for general troubleshooting. + +## References +Hersbach, H., Bell, B., Berrisford, P., Biavati, G., Horányi, A., Muñoz Sabater, J., Nicolas, J., Peubey, C., Radu, R., Rozum, I., Schepers, D., Simmons, A., Soci, C., Dee, D., Thépaut, J-N. (2023): ERA5 hourly data on pressure levels from 1940 to present. Copernicus Climate Change Service (C3S) Climate Data Store (CDS), DOI: 10.24381/cds.bd0915c6 (Accessed on 11 Jan. 2026) \ No newline at end of file diff --git a/data_access/homework2/control_flow.py b/data_access/homework2/control_flow.py new file mode 100644 index 0000000..34ada9f --- /dev/null +++ b/data_access/homework2/control_flow.py @@ -0,0 +1,78 @@ +import argparse +import random +from datetime import datetime, timedelta, date + +# Define start and end date for download routine +START_DATE = date(2026, 1, 1) # Begin of project +ARCHIVE = set() # mocks a folder of finished days + +def download_data(day): + """ + This function "downloads" data and fails randomly. + """ + if random.random() < 0.3: # 30% failure rate + raise ConnectionError(f"Download failed for {day}") + print(f" [1/3] Downloaded data for {day}") + +def process_data(day): + """ + Process test data. + """ + print(f" [2/3] Processed data for {day}") + +def archive_data(day): + """ + Add downloaded data to archive. + """ + ARCHIVE.add(day) + print(f" [3/3] Archived {day}") + +def run_pipeline(day): + """ + This function contains the whole control flow for downloading data. + In this case nothing is really downloaded. + """ + print(f"\n--- Starting Pipeline for {day} ---") + try: + if day in ARCHIVE: + print(f" Skipping {day}: Already archived.") + return True + + download_data(day) + process_data(day) + archive_data(day) + return True + except Exception as e: + print(f" ERROR: {e}") + return False + +# need to call pipeline function +def main(): + parser = argparse.ArgumentParser(description="Daily Data Manager") + parser.add_argument("--date", type=lambda s: datetime.strptime(s, "%Y-%m-%d").date(), + help="Process a specific date (YYYY-MM-DD)") + args = parser.parse_args() + + if args.date: + # Scenario 1: Specific date requested + run_pipeline(args.date) + else: + # Scenario 2: Find gaps and process until today + today = date.today() + current = START_DATE + print(f"Auto-syncing from {START_DATE} to {today}...") + + while current <= today: + if current not in ARCHIVE: + success = run_pipeline(current) + if not success: + print(f"Stopping at {current} due to failure.") + break + current += timedelta(days=1) + +if __name__ == "__main__": + # Pre-filling archive with some "finished" dates to test gap detection + ARCHIVE.add(date(2026, 1, 1)) + ARCHIVE.add(date(2026, 1, 2)) + + main() \ No newline at end of file diff --git a/data_access/homework2/load_ERA5_data.py b/data_access/homework2/load_ERA5_data.py new file mode 100644 index 0000000..f53c2e4 --- /dev/null +++ b/data_access/homework2/load_ERA5_data.py @@ -0,0 +1,17 @@ +### IMPORT MODULES +import ERA5_control_flow as cf + + +### PARAMETER TO CHANGE ### + +VARIABLE = ["specific_humidity"] # change variable here +TIME = ["00:00", "06:00", "12:00","18:00"] # download 6-hourly data +PRESSURE_LEVEL = ["300", "500", "800", "900", "975"] # change pressure levels here + +start = "2024-12-01" +end = "2024-12-05" + +DATASET = "reanalysis-era5-pressure-levels" + + +cf.retrieve_and_process(DATASET, start, end, VARIABLE, TIME, PRESSURE_LEVEL) diff --git a/data_access/homework2/plot_ERA5_data.py b/data_access/homework2/plot_ERA5_data.py new file mode 100644 index 0000000..b4712a4 --- /dev/null +++ b/data_access/homework2/plot_ERA5_data.py @@ -0,0 +1,93 @@ +import xarray as xr +import healpy as hp +import matplotlib.pyplot as plt +import numpy as np +import cartopy.crs as ccrs +import cartopy.feature as cfeature + +ds_hp = xr.open_zarr("era5_nside16.zarr", consolidated=True) + +times = ds_hp.time.values[[2, 7]] +sample = ds_hp.sel(time=times) + +## plot 1: regridded data +data = sample['q'].sel(time=times[0], level=500).values +hp.mollview(data, title=f"Specific humidity 500 hPa\n{times[0]}") +plt.show() + +## plot 2: regridded data +data = sample['q'].sel(time=times[1], level=500).values +hp.mollview(data, title=f"Specific humidity 500 hPa\n{times[1]}") +plt.show() + +#### Plot original netCDF file +ds_nc1 = xr.open_dataset('era5_2024-12-01.nc', engine='netcdf4') +humidity1 = ds_nc1['q'] + +ds_nc2 = xr.open_dataset('era5_2024-12-02.nc', engine='netcdf4') +humidity2 = ds_nc2['q'] + +# Select data +t_data1 = humidity1.sel( + pressure_level=500, + valid_time=times[0] +) + +t_data2 = humidity2.sel( + pressure_level=500, + valid_time=times[1] +) + +# Create the plot with Robinson projection +fig, (ax1, ax2) = plt.subplots(nrows=1,ncols=2,figsize=(16, 4),subplot_kw={'projection': ccrs.Robinson()}) + +# Plot the temperature field +cf1 = t_data1.plot( + ax=ax1, + transform=ccrs.PlateCarree(), # Data is in PlateCarree (lon/lat) + cmap='coolwarm', + cbar_kwargs={'label': 'humidity (kg/kg)'}, + add_colorbar=True, + add_labels=False # Avoid duplicate labels +) + +# Set title +ax1.set_title('Plot with netCDF data', fontsize=14, pad=20) + +# Add coastlines and country borders +ax1.add_feature(cfeature.COASTLINE, linewidth=0.8) +ax1.add_feature(cfeature.BORDERS, linewidth=0.5, linestyle=':') +ax1.add_feature(cfeature.LAND, color='lightgray') +ax1.add_feature(cfeature.OCEAN, color='lightblue') + +# Set global view +ax1.set_global() + +### for 2nd plot +# Plot the temperature field +cf2 = t_data2.plot( + ax=ax2, + transform=ccrs.PlateCarree(), # Data is in PlateCarree (lon/lat) + cmap='coolwarm', + cbar_kwargs={'label': 'humidity (kg/kg)'}, + add_colorbar=True, + add_labels=False # Avoid duplicate labels +) + +# Set title +ax2.set_title('Plot with netCDF data', fontsize=14, pad=20) + +# Add coastlines and country borders +ax2.add_feature(cfeature.COASTLINE, linewidth=0.8) +ax2.add_feature(cfeature.BORDERS, linewidth=0.5, linestyle=':') +ax2.add_feature(cfeature.LAND, color='lightgray') +ax2.add_feature(cfeature.OCEAN, color='lightblue') + +# Set global view +ax2.set_global() + + +# Improve layout +plt.tight_layout() + +plt.show() \ No newline at end of file diff --git a/data_access/load_noaa_class_sst.ipynb b/data_access/load_noaa_class_sst.ipynb new file mode 100644 index 0000000..28e8bc0 --- /dev/null +++ b/data_access/load_noaa_class_sst.ipynb @@ -0,0 +1,1363 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "5a03c248-1f5f-485e-801a-ca121cd5d427", + "metadata": {}, + "source": [ + "# Homework 1 – Data Access Notebook\n", + "**Author:** Johanna Kasischke \n", + "**Course:** ESDP \n", + "**Due Date:** 04.12.2025" + ] + }, + { + "cell_type": "markdown", + "id": "da5ff78a-5512-4a2d-aef3-0e049bf11865", + "metadata": {}, + "source": [ + "This notebook demonstrates how to access and download NOAA NESDIS SST14NA data, estimate the data volume, and open the NetCDF file for analysis." + ] + }, + { + "cell_type": "markdown", + "id": "2a929d1a-226a-4df9-89d0-5c685f84a2f0", + "metadata": {}, + "source": [ + "## NESDIS 14 km sea surface temperature for North America, 1993 to 2016" + ] + }, + { + "cell_type": "markdown", + "id": "892b609b-a5ba-4b9e-9b4e-5f44b57a5874", + "metadata": {}, + "source": [ + "CLASS (Comprehensive Large Array-data Stewardship System) is an electronic library of NOAA environmental data. This website provides the data from the POES (Polar Orbiting Operational Environmental Satellite) and GOES (Geostationary Operational Environmental Satellite) satellite programmes (https://ospo.noaa.gov/resources/data-access/).\n", + "The National Oceanic and Atmospheric Administration (NOAA) is the United States' federal agency for oceanic and atmospheric research. The dataset used was distributed by the National Environmental Satellite, Data and Information Service (NESDIS). This organisation operates and manages US satellite services. One of its tasks is to monitor environmental conditions on a daily, weekly, seasonal and yearly basis. " + ] + }, + { + "cell_type": "markdown", + "id": "f81af9f9-7898-4007-9d73-11c1cf76154f", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "## Introduction to dataset\n", + "\n", + "This task uses a dataset titled 'NESDIS 14 km sea surface temperature for North America, 1993 to 2016'. It shows sea surface temperatures in degrees Celsius. Sea surface temperature (SST) is the temperature of the ocean's surface water. The time period of interest is the Atlantic hurricane season, from 8 June 2025 to 6 January 2006, because SST has a significant impact on hurricane development." + ] + }, + { + "cell_type": "markdown", + "id": "04d8ac02-55b8-415e-9788-db3bc5fe5ded", + "metadata": {}, + "source": [ + "## How to download the data" + ] + }, + { + "cell_type": "markdown", + "id": "90f1820e-fc64-4f12-8186-87396f6a2101", + "metadata": {}, + "source": [ + "The data can be downloaded from the NOAA Comprehensive Large Array-data Stewardship System (CLASS), which is an electronic library of NOAA environmental data.\n", + "\n", + "The NESDIS SST14NA dataset was accessed via: https://www.aev.class.noaa.gov/saa/products/welcome. The following steps need to be followed to download the data:\n", + "1. Navigate to the portal and select \"Sea surface temperature (14 km North America)\" either at the top search bar or on the right-hand side of the website.\n", + "2. Specify the temporal range (here: June 8 2005 till January 6 2006).\n", + "3. Choose output format (NetCDF preferred). Other options include ASCII, GIF, spreadsheet, tab-delimited or CSV.\n", + "4. Select the distributed variable (SST).\n", + "5. Initiate the search --> the server retrieves the specific data\n", + "6. The download starts by itself." + ] + }, + { + "cell_type": "markdown", + "id": "d999aa8e-21cc-49c6-9d20-2ac8af56c392", + "metadata": {}, + "source": [ + "*Citation of the dataset:*\n", + "NOAA/NESDIS Office of Satellite Data Processing and Distribution (2013). NESDIS 14 km sea surface temperature for North America, 1993 to 2016. Data used from 2005-06-08 until 2006-01-06. NOAA National Centers for Environmental Information. Dataset. https://www.ncei.noaa.gov/archive/accession/OSDPD-SST-14km-NorthAmerica. Accessed 30th of November 2025." + ] + }, + { + "cell_type": "markdown", + "id": "c996d2fb-14dc-485b-9f50-c3102b12de6e", + "metadata": {}, + "source": [ + "### Estimate of data volume for download" + ] + }, + { + "cell_type": "markdown", + "id": "558e9ec8-8cbc-4013-b29b-48d675f8ba1b", + "metadata": {}, + "source": [ + "- spatial domain: North America ((-140° W, -50° E, 15° S, 60° N)\n", + "- Spatial resolution: 14 km\n", + "- Pixel resolution: .125 x .125 degrees\n", + "- Time range: 8th of June, 2025 until 6th of January, 2006 (Atlantic Hurricane Season of 2005)\n", + "- Time resolution: 48 h" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "d57dd8e9-4460-4600-ad9a-407c9c45993e", + "metadata": {}, + "outputs": [], + "source": [ + "import datetime\n", + "import math\n", + "import numpy as np\n", + "\n", + "\n", + "# Here, you can change the parameters for specific data set and then compute the approximate download volume.\n", + "\n", + "spatial_resolution = 14 # [km]\n", + "degree_resolution = 0.125 # The pixel resolution is .125 x .125 degrees. \n", + "lat_min, lat_max = -15, 60\n", + "lon_min, lon_max = -140, -50\n", + "date_start = datetime.date(2005, 6, 8)\n", + "date_end = datetime.date(2006, 1, 6)\n", + "time_step = 2 # equals the 48 h time resolution\n", + "bytes = 4 # assume 4 bytes/value" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "c0fc4672-c259-40b0-a454-491c110c5708", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "For calculation with spatial resolution, the estimated download volume will be: 169.0 mb\n", + "For calculation with pixel resolution, the estimated download volume will be: 185.0 mb\n" + ] + } + ], + "source": [ + "days = (date_end - date_start).days\n", + "steps = days / time_step + 1\n", + "\n", + "degree_to_km = 111.32 # approx. km per degree latitude\n", + "\n", + "cell_height = (lat_max - lat_min) * degree_to_km\n", + "\n", + "# for approximating the download size, the latitude will be approximated for the whole calculation (-> longitude distance varies with latitude)\n", + "lat_mean = (lat_min + lat_max) / 2\n", + "lon_per_degree = degree_to_km * math.cos(math.radians(lat_mean)) # in km\n", + "cell_width = (lon_max - lon_min) * lon_per_degree\n", + "\n", + "# calculate grid size\n", + "nx = (cell_width / spatial_resolution)\n", + "ny = (cell_height / spatial_resolution)\n", + "pixels = nx * ny\n", + "\n", + "# with pixel resolution\n", + "nx_p = (lon_max - lon_min) / degree_resolution\n", + "ny_p = (lat_max - lat_min) / degree_resolution\n", + "pixels_p = nx_p * ny_p\n", + "\n", + "\n", + "# compute total bytes with spatial resolution\n", + "total_bytes = pixels * bytes * steps\n", + "total_mb = total_bytes / 1e6\n", + "\n", + "\n", + "# compute total bytes with pixel resolution\n", + "total_bytes_p = pixels_p * bytes * steps\n", + "total_mb_p = total_bytes_p / 1e6\n", + "\n", + "\n", + "########## PRINT SOLUTION #############\n", + "print(\"For calculation with spatial resolution, the estimated download volume will be: \", np.round(total_mb), \"mb\")\n", + "print(\"For calculation with pixel resolution, the estimated download volume will be: \", np.round(total_mb_p), \"mb\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "e0fccc64-5acc-4fe3-9785-ba803e8ab373", + "metadata": {}, + "source": [ + "Note that the estimated and real data volumes don't match. This might be due to the data stored as netCDF. " + ] + }, + { + "cell_type": "markdown", + "id": "4601e454-77f9-4e9b-84e4-4d198301b53b", + "metadata": {}, + "source": [ + "## Download and import the data" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "a554dac9-f2a1-4ca9-a58d-80da827551b5", + "metadata": {}, + "outputs": [], + "source": [ + "# import modules\n", + "import xarray as xr" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "id": "70183c02-92a5-45bf-a54a-a2d0926d0ae5", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Dimensions: (AX005: 226, LAT: 361, TIME1: 107, bnds: 2)\n", + "Coordinates:\n", + " * AX005 (AX005) float64 220.0 220.4 220.8 221.2 ... 309.2 309.6 310.0\n", + " * LAT (LAT) float32 15.0 15.12 15.25 15.38 ... 59.62 59.75 59.88 60.0\n", + " * TIME1 (TIME1) datetime64[ns] 2005-06-09 2005-06-11 ... 2006-01-06\n", + "Dimensions without coordinates: bnds\n", + "Data variables:\n", + " TIME1_bnds (TIME1, bnds) datetime64[ns] ...\n", + " ANAL_TEMP (TIME1, LAT, AX005) float32 ...\n", + "Attributes:\n", + " history: FERRET V6.96 30-Nov-25\n", + " Conventions: CF-1.6\n", + "\n", + "Dataset opened successfully\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/srv/conda3/envs/py39-cartopy/lib/python3.9/site-packages/xarray/conventions.py:431: SerializationWarning: variable 'ANAL_TEMP' has multiple fill values {9999.0, 999.9}, decoding all values to NaN.\n", + " new_vars[k] = decode_cf_variable(\n", + "/srv/conda3/envs/py39-cartopy/lib/python3.9/site-packages/xarray/conventions.py:431: SerializationWarning: variable 'ANAL_TEMP' has multiple fill values {9999.0, 999.9}, decoding all values to NaN.\n", + " new_vars[k] = decode_cf_variable(\n" + ] + } + ], + "source": [ + "# upload file into ipynb\n", + "files = [\"2005-06-08_00_00_00_2006-01-06_23_59_59-SST14NA-anal_temp-15_60_220_310-.5.NETCDF\"]\n", + "\n", + "# open the file in xarray\n", + "datasets = []\n", + "for file in files:\n", + " ds = xr.open_dataset(file)\n", + " print(ds)\n", + " datasets.append(ds)\n", + "try:\n", + " ds = xr.open_dataset(file)\n", + " print(\"\\nDataset opened successfully\")\n", + "except Exception as e:\n", + " print(\"Error opening dataset:\", e)\n" + ] + }, + { + "cell_type": "markdown", + "id": "84cc3083-e7d0-41b0-b0e0-bfaa6ccbbdd9", + "metadata": {}, + "source": [ + "#### Verify download\n", + "In order to verify that each dataset has been successfully downloaded, the length of one array and the expected time steps will be compared. The estimation of the download volume was achieved through the computation of time steps. The following Python code verifies that each dataset has been downloaded. Therefore the length of the datetime array and the number of time steps must be equal." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "06923225-78c8-4639-8f00-fc87e2b867d7", + "metadata": {}, + "outputs": [], + "source": [ + "len(ds[\"TIME1_bnds\"].values)" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "0b2b3737-ee1e-4311-a2d7-8238c1f8b5c1", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Frozen({'AX005': \n", + "array([220. , 220.4, 220.8, ..., 309.2, 309.6, 310. ])\n", + "Attributes:\n", + " units: degrees_east\n", + " point_spacing: even\n", + " axis: X\n", + " standard_name: longitude, 'LAT': \n", + "array([15. , 15.125, 15.25 , ..., 59.75 , 59.875, 60. ], dtype=float32)\n", + "Attributes:\n", + " long_name: Latitude\n", + " units: degrees_north\n", + " point_spacing: even\n", + " axis: Y\n", + " standard_name: latitude, 'TIME1': \n", + "array(['2005-06-09T00:00:00.000000000', '2005-06-11T00:00:00.000000000',\n", + " '2005-06-13T00:00:00.000000000', '2005-06-15T00:00:00.000000000',\n", + " '2005-06-17T00:00:00.000000000', '2005-06-19T00:00:00.000000000',\n", + " '2005-06-21T00:00:00.000000000', '2005-06-23T00:00:00.000000000',\n", + " '2005-06-25T00:00:00.000000000', '2005-06-27T00:00:00.000000000',\n", + " '2005-06-29T00:00:00.000000000', '2005-07-01T00:00:00.000000000',\n", + " '2005-07-03T00:00:00.000000000', '2005-07-05T00:00:00.000000000',\n", + " '2005-07-07T00:00:00.000000000', '2005-07-09T00:00:00.000000000',\n", + " '2005-07-11T00:00:00.000000000', '2005-07-13T00:00:00.000000000',\n", + " '2005-07-15T00:00:00.000000000', '2005-07-17T00:00:00.000000000',\n", + " '2005-07-19T00:00:00.000000000', '2005-07-21T00:00:00.000000000',\n", + " '2005-07-23T00:00:00.000000000', '2005-07-25T00:00:00.000000000',\n", + " '2005-07-27T00:00:00.000000000', '2005-07-29T00:00:00.000000000',\n", + " '2005-07-31T00:00:00.000000000', '2005-08-02T00:00:00.000000000',\n", + " '2005-08-04T00:00:00.000000000', '2005-08-06T00:00:00.000000000',\n", + " '2005-08-08T00:00:00.000000000', '2005-08-10T00:00:00.000000000',\n", + " '2005-08-12T00:00:00.000000000', '2005-08-14T00:00:00.000000000',\n", + " '2005-08-16T00:00:00.000000000', '2005-08-18T00:00:00.000000000',\n", + " '2005-08-20T00:00:00.000000000', '2005-08-22T00:00:00.000000000',\n", + " '2005-08-24T00:00:00.000000000', '2005-08-26T00:00:00.000000000',\n", + " '2005-08-28T00:00:00.000000000', '2005-08-30T00:00:00.000000000',\n", + " '2005-09-01T00:00:00.000000000', '2005-09-03T00:00:00.000000000',\n", + " '2005-09-05T00:00:00.000000000', '2005-09-07T00:00:00.000000000',\n", + " '2005-09-09T00:00:00.000000000', '2005-09-11T00:00:00.000000000',\n", + " '2005-09-13T00:00:00.000000000', '2005-09-15T00:00:00.000000000',\n", + " '2005-09-17T00:00:00.000000000', '2005-09-19T00:00:00.000000000',\n", + " '2005-09-21T00:00:00.000000000', '2005-09-23T00:00:00.000000000',\n", + " '2005-09-25T00:00:00.000000000', '2005-09-27T00:00:00.000000000',\n", + " '2005-09-29T00:00:00.000000000', '2005-10-01T00:00:00.000000000',\n", + " '2005-10-03T00:00:00.000000000', '2005-10-05T00:00:00.000000000',\n", + " '2005-10-07T00:00:00.000000000', '2005-10-09T00:00:00.000000000',\n", + " '2005-10-11T00:00:00.000000000', '2005-10-13T00:00:00.000000000',\n", + " '2005-10-15T00:00:00.000000000', '2005-10-17T00:00:00.000000000',\n", + " '2005-10-19T00:00:00.000000000', '2005-10-21T00:00:00.000000000',\n", + " '2005-10-23T00:00:00.000000000', '2005-10-25T00:00:00.000000000',\n", + " '2005-10-27T00:00:00.000000000', '2005-10-29T00:00:00.000000000',\n", + " '2005-10-31T00:00:00.000000000', '2005-11-02T00:00:00.000000000',\n", + " '2005-11-04T00:00:00.000000000', '2005-11-06T00:00:00.000000000',\n", + " '2005-11-08T00:00:00.000000000', '2005-11-10T00:00:00.000000000',\n", + " '2005-11-12T00:00:00.000000000', '2005-11-14T00:00:00.000000000',\n", + " '2005-11-16T00:00:00.000000000', '2005-11-18T00:00:00.000000000',\n", + " '2005-11-20T00:00:00.000000000', '2005-11-22T00:00:00.000000000',\n", + " '2005-11-24T00:00:00.000000000', '2005-11-26T00:00:00.000000000',\n", + " '2005-11-28T00:00:00.000000000', '2005-11-30T00:00:00.000000000',\n", + " '2005-12-02T00:00:00.000000000', '2005-12-04T00:00:00.000000000',\n", + " '2005-12-06T00:00:00.000000000', '2005-12-08T00:00:00.000000000',\n", + " '2005-12-10T00:00:00.000000000', '2005-12-12T00:00:00.000000000',\n", + " '2005-12-14T00:00:00.000000000', '2005-12-16T00:00:00.000000000',\n", + " '2005-12-18T00:00:00.000000000', '2005-12-20T00:00:00.000000000',\n", + " '2005-12-22T00:00:00.000000000', '2005-12-24T00:00:00.000000000',\n", + " '2005-12-26T00:00:00.000000000', '2005-12-28T00:00:00.000000000',\n", + " '2005-12-30T00:00:00.000000000', '2005-12-31T00:00:00.000000000',\n", + " '2006-01-02T00:00:00.000000000', '2006-01-04T00:00:00.000000000',\n", + " '2006-01-06T00:00:00.000000000'], dtype='datetime64[ns]')\n", + "Attributes:\n", + " time_origin: 1-JAN-1997 00:00:00\n", + " axis: T\n", + " standard_name: time\n", + " bounds: TIME1_bnds, 'TIME1_bnds': \n", + "[214 values with dtype=datetime64[ns]], 'ANAL_TEMP': \n", + "[8729702 values with dtype=float32]\n", + "Attributes:\n", + " long_name: Analysis Temperature\n", + " units: Deg. C\n", + " long_name_mod: regrid: 0.4 deg on X\n", + " history: From CLASS_Descriptor})\n", + "Frozen({'AX005': 226, 'LAT': 361, 'TIME1': 107, 'bnds': 2})\n", + "\n", + "[8729702 values with dtype=float32]\n", + "Coordinates:\n", + " * AX005 (AX005) float64 220.0 220.4 220.8 221.2 ... 308.8 309.2 309.6 310.0\n", + " * LAT (LAT) float32 15.0 15.12 15.25 15.38 ... 59.62 59.75 59.88 60.0\n", + " * TIME1 (TIME1) datetime64[ns] 2005-06-09 2005-06-11 ... 2006-01-06\n", + "Attributes:\n", + " long_name: Analysis Temperature\n", + " units: Deg. C\n", + " long_name_mod: regrid: 0.4 deg on X\n", + " history: From CLASS_Descriptor\n" + ] + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 20, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj4AAAHFCAYAAADyj/PrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAACd2UlEQVR4nO3dd3gUVfcH8O+kbQpJIIQ0CAkdQjdBCCAJXQRpKiqoRJAXpIYuihAQQu/+AFFpKuqrgPKKQKI0EdAQQEKRopQAiREIaaTv/f0RZtjd7G5mZvvu+TzPPpDZ2bt3J2XPnnvuvRxjjIEQQgghxAE4WboDhBBCCCHmQoEPIYQQQhwGBT6EEEIIcRgU+BBCCCHEYVDgQwghhBCHQYEPIYQQQhwGBT6EEEIIcRgU+BBCCCHEYVDgQwghhBCHQYEPsXocx4m6HT58GDdu3ADHcVi+fLnw+MOHDwvnbN26VetzdOvWDRzHITw8XO14eHi4zueLjY0VzsvLy8OMGTPQq1cv1KpVCxzHISEhwfgXQ4SDBw9ixIgRaNq0Kby8vFC7dm0MGDAAqampWs8/ffo0evTogWrVqqF69eoYPHgw/v77b63nrlu3Dk2bNoVCoUC9evUwb948lJaWqp2zdetWndcsMzNT9OuQ0q+bN29ixIgRCAkJgUKhQO3atTFo0CDRz5Wfn4/4+HiEhITA3d0dbdq0wVdffVXpPMYY1q5dK1yD4OBgvP3228jOzhb9XIQQy3KxdAcIqcqJEyfUvv7ggw9w6NAhHDx4UO14REQEHjx4oLMdb29vfPrpp4iLi1M7fv36dRw+fBg+Pj5aH9epUye1QIqnev79+/exadMmtG7dGgMHDsQnn3xS1csymQ0bNuD+/fuYNGkSIiIi8O+//2LFihXo0KEDDhw4gG7dugnn/vnnn4iNjUWbNm3w3//+F0VFRZgzZw6eeeYZnD17FrVq1RLOXbhwId5//32888476NWrF1JSUjB79mzcuXMHmzZtqtSPLVu2oGnTpmrHatasKeo1SOnX+fPnERsbi/r162P58uWoU6cOMjIycODAAdHXbPDgwUhJScHixYvRuHFj7NixA6+++iqUSiWGDh0qnDdt2jSsXr0a06ZNQ48ePXDx4kXMmTMHKSkpOHHiBFxdXUU/JyHEQhghNmb48OHMy8tL633Xr19nANiyZcuEY4cOHWIA2FtvvcUAsCtXrqg9Zvbs2axOnTqsT58+LCwsTO2+sLAw1rdv3yr7pFQqmVKpZIwx9u+//zIAbO7cudJemJH8888/lY7l5eWxwMBA1r17d7XjL730EvP392c5OTnCsRs3bjBXV1c2Y8YM4di9e/eYu7s7+89//qP2+IULFzKO49iFCxeEY1u2bGEAWEpKiuzXILZfSqWStWnThrVp04YVFRXJeq69e/cyAGzHjh1qx3v27MlCQkJYWVkZY4yx27dvM2dnZzZhwgS183bs2MEAsE2bNsl6fkKIedFQF3EYPXv2RGhoKDZv3iwcUyqV2LZtG4YPHw4nJ/m/DvxQjjUICAiodKxatWqIiIhAenq6cKysrAw//PADXnjhBbXsVVhYGLp27Yrdu3cLx/bv34+ioiK8+eabau2++eabYIzhu+++M1r/pfTr6NGjOHv2LOLj46FQKGQ93+7du1GtWjW89NJLasfffPNN3L17F7/99hsA4OTJkygvL8dzzz2ndl6/fv0AADt37pT1/IQQ86LAhzgMJycnxMXFYfv27SgvLwcAJCUl4fbt25Xe0FUxxlBWVlbpxhgzWt90PYe2mxw5OTk4ffo0mjdvLhz766+/UFhYiFatWlU6v1WrVrh27RqKiooAVAwnAUDLli3VzgsODoa/v79wv6p+/frB2dkZfn5+GDx4sNZztJHSr6NHjwKoGMZ87rnn4O7ujmrVqqFfv374888/RT3f+fPn0axZM7i4qI/888/P97ukpAQAKgVYrq6u4DgO586dE/V8hBDLosCHOJQ333wTGRkZ2L9/PwBg8+bNiImJQYMGDXQ+5scff4Srq2ul28KFC43Wr23btml9Dm03OcaNG4eCggK89957wrH79+8DAPz8/Cqd7+fnB8aYULR7//59KBQKeHl5aT2XbwsAgoKC8N577+GTTz7BoUOH8MEHHyAlJQUdOnTAH3/8UWVfpfTrzp07ACq+ryEhIdi7dy82btyI8+fP45lnnkFGRoao59P1XKr9iYiIAAD8+uuvaucdP34cjDG1a0AIsV5U3EwcSr169RAbG4vNmzejQ4cO+P7776ssRO7cuTNWrVpV6Xjt2rWN1q/nn38eKSkpRmtP1fvvv48vvvgC69atQ2RkZKX79Q3Rqd4n9rxnn30Wzz77rPB1ly5d0LdvX7Rs2RJz5szB999/D6BimFGpVKq14ezsLOn5+MdHR0erfR9btGiBtm3b4v/+7/+wYMECMMaELB9PNcMj5rlat26NLl26YNmyZWjSpAl69uyJixcvYsyYMXB2djZoqJQQYj4U+BCHM3LkSLz55ptYuXIlPDw88OKLL+o939fXF1FRUSbtk5+fH3x9fY3e7rx587BgwQIsXLgQ48ePV7uPn2GlLVPx4MEDcByH6tWrC+cWFRXh0aNH8PT0rHSutoBKVXh4ODp37oyTJ08Kx0aMGIFt27YJX8fExODw4cOS+wUAvXv3VjuvTZs2CA4OxunTpwEAR44cQdeuXdXOuX79OsLDw1GzZk2dzwWoZ56++eYbxMXFYciQIQAANzc3TJ48GT/99BMePnyo9xoQQqwDBT7E4QwePBjjxo3D4sWLMWrUKHh4eFi6S9i2bZveOiNVYmuL5s2bh4SEBCQkJODdd9+tdH+DBg3g4eGBtLS0SvelpaWhYcOGcHd3B/CktictLQ3t27cXzsvMzMS9e/fQokULUf1WzYokJCSoBWPe3t6S+6WtDkjb80VGRlbKqIWEhAiv7csvv0RZWZlaFoh/ftXXFhAQgB9//BFZWVnIzMxEWFgYPDw8sH79+ioDaEKIdaDAhzgcDw8PzJkzB0ePHsXbb79t6e4AMP5Q1wcffICEhATMnj0bc+fO1XqOi4sLnn/+eezatQtLly4VAo9bt27h0KFDmDx5snDus88+C3d3d2zdulUt8OEXKxw4cKDe/ly/fh2//vorevToIRwLDw+vtGCk1H716dMHnp6e2Ldvn9rx06dPIzMzEx06dABQEVTpytoNGjQIH3/8MXbu3ImXX35ZOL5t2zaEhISovV5eQECAMHtu7dq1KCgoqJRRI4RYJwp8iEOaMmUKpkyZIurchw8fqg3R8BQKBdq2bSt8vW/fPhQUFCAvLw8AcPHiRXz77bcAgOeee67SEJGqmjVril7cryorVqzAnDlz8Oyzz6Jv376V+s4HA0BFVqhdu3bo168f3nnnHWGhQH9/f0ydOlU4z8/PD7Nnz8b7778PPz8/YQHDhIQEvPXWW0LhLwD06NEDXbp0QatWreDj44O0tDQsXboUHMfhgw8+EPUaxParevXqmD9/PqZNm4a4uDi8+uqryMzMxPvvv4+6deti7NixVT5Xnz590LNnT7z99tvIzc1Fw4YN8eWXX2L//v34/PPP1eqOPv74YwAVWamHDx9i3759+PTTT5GYmIinnnpK1GsjhFiYhdYPIkQ2uQsYfvPNN3rb7du3r9YFDAFovdWuXVv0udevX5f1WuWIiYnR2Q9tv/KnTp1i3bt3Z56enszHx4cNHDiQXbt2TWvba9asYY0bN2Zubm6sbt26bO7cuaykpETtnPj4eBYREcG8vb2Zi4sLCwkJYa+99hq7fPmypNchpV8ff/wxa9GiBXNzc2M1a9Zkw4YNY+np6aKfKy8vj02cOJEFBQUxNzc31qpVK/bll19WOu+jjz5izZo1Y56enqxatWrsmWeeYd99952k10UIsSyOMSMuRkIIIYQQYsVo/iUhhBBCHAYFPoQQQghxGBT4EEIIIcRhWDzwuXPnDl577TXUrFkTnp6eaNOmDVJTU4X7GWNISEhASEgIPDw8EBsbiwsXLliwx4QQQgixVRYNfLKzs9GpUye4urpi3759uHjxIlasWCGsygoAS5cuxcqVK/Hhhx8iJSUFQUFB6NmzpzBlmBBCCCFVKysrw+zZs1GvXj14eHigfv36mD9/vtrWMWKSDcXFxZgwYQL8/f3h5eWF/v374/bt2+Z+OfJZckrZzJkzWefOnXXer1QqWVBQEFu8eLFwrKioiPn6+rKNGzeao4uEEEKIXViwYAGrWbMm++GHH9j169fZN998w6pVq8ZWr14tnLN48WLm7e3Ndu7cydLS0tjLL7/MgoODWW5urnDOmDFjWO3atVlycjI7ffo069q1K2vdujUrKyuzxMuSzKLT2SMiItC7d2/cvn0bR44cQe3atTF27FiMGjUKAPD333+jQYMGOH36tNpCcQMGDED16tXV9vnhFRcXo7i4WPhaqVTiwYMHqFmzpt6NCAkhhBDGGPLy8hASEmLSjWeLiopQUlJicDtubm7CFi5V6devHwIDA/Hpp58Kx1544QV4enris88+A2MMISEhiI+Px8yZMwFUvKcGBgZiyZIlGD16NHJyclCrVi189tlnwkrnd+/eRWhoKH788cdK++ZZJUtGXQqFgikUCjZr1ix2+vRptnHjRubu7s62bdvGGGPs119/ZQDYnTt31B43atQo1qtXL61tzp07V+/ibXSjG93oRje6VXWTsgCmVIWFhcwTzkbpZ1BQEPvnn39YTk6OcCsqKtL6vIsWLWJhYWHCYqJnz55lAQEBbMeOHYwxxv766y8GgJ0+fVrtcf3792dvvPEGY4yxn3/+mQFgDx48UDunVatWbM6cOca+VCZh0S0rlEoloqKikJiYCABo27YtLly4gA0bNuCNN94QztPM1DDGdGZvZs2apbYVQU5ODurWrYv09HTMCe0IAFidc97YL4UQm1VUWGjpLliUu4cHCouKdN7PaSTF3a1gU1uevn6LwTFW5etxtJ+PvLw8NGzUSNgjzhRKSkrwCOV4A7XhZkCpbQmU2J55B4GBgWrH586di4SEhErnz5w5Ezk5OWjatCmcnZ1RXl6OhQsX4tVXXwVQsekwgErtBQYG4ubNm8I5bm5uqFGjRqVz+MdbO4sGPsHBwWp7/ABAs2bNsHPnTgBAUFAQgIoLHRwcLJyTlZVV6RvDUygUUCgUlY77+PhgM7tlrK4TYjfcXF0t3QWLcvfwEK4B0/KBig98rCng4bm6uVU6xvdX22vR5CFiiETbzwfftmZQaE/MURrhBie4cQYMpz2+/Onp6fDx8REOa3sPBICvv/4an3/+OXbs2IHmzZvj7NmziI+PR0hICIYPHy6cJyXZIOUca2HRWV2dOnXC5cuX1Y5duXIFYWFhAIB69eohKCgIycnJwv0lJSU4cuQIOnbsKOm5ih4VGN5hQojdiPdshnjPZmoZDY4xtTdza39j5/uretO8T5+yM/sN7gPjOFFBFqnMmeMMvgEVH+xVb7oCn+nTp+Odd97BK6+8gpYtW+L111/H5MmTsWjRIgDqyQZVqsmGoKAglJSUIDs7W+c51s6igc/kyZNx8uRJJCYm4tq1a9ixYwc2bdqEcePGAaiIOuPj45GYmIjdu3fj/PnziIuLg6enJ4YOHWrJrhM7MIYLx90Fb+Pugrct3RViAR+e3ogPT2+0dDdkGcOFYwwXbnA7Lm2flfU41aBKTIBFtHPiAGcDbk4S481Hjx5VKth2dnYWprOLSTZERkbC1dVV7ZyMjAycP39eckLCUiw61NWuXTvs3r0bs2bNwvz581GvXj2sXr0aw4YNE86ZMWMGCgsLMXbsWGRnZ6N9+/ZISkqSNf5amrIHAODarr/RXgOxXasfXap0TPXTvzUObRDzsPY38o3sBgD99TfGysK4e3g4XJ2PvXr++eexcOFC1K1bF82bN8eZM2ewcuVKjBgxAoB6sqFRo0Zo1KgREhMT1ZINvr6+GDlyJKZOnYqaNWvCz88P06ZNQ8uWLdGjRw9LvjzR7H539tzcXPj6+iLr7m1hDPTe8mmoPfcjC/eMWFpVf8wdJfAx5ZuaPdSCWPPPwU/N2wEAOp86Wula6wt8pNYtOVLgk5ubi8CgIOTk5KjVzRj7OXx9fTHeOQwKA2p8ipkSH5bfFN3XvLw8vP/++9i9ezeysrIQEhKCV199FXPmzIHb43oxxhjmzZuHjz76SEg2/N///R9atGghtFNUVITp06djx44dKCwsRPfu3bF+/XqEhobKfi3m5JCBjyqFl+kq94l1oqyOOs03NTFvmI7EFn5GtAUmxipuBtQDLJvElIDIAMOcgc8kF8MDnzVl4gMfUsGiQ13WoLigYusLCoDsH//mYAtvZMQ62PLPirYgVe7wV48LKQBsOPNjyMwpYnccPvAhjiPesxmAJ/URpILLn0cAAKXNYi3bESKbu4eHMEOrrGmMzvMmezQFQL8D1kJ1Zpasx4Nm08nhMIEPqyLip8yP/eI/pdIfe+34mT3a/hgYukAeMR/++1imJytj6O/AgxUVi8P6TV1pUDukAj87S/bjjdcVh+IwgQ/wJPjhmLKKM7UrelQAd08vY3aJGBm/XpPq98mWhyuIZRyL6gLgyRAPIcR+OFTgw2Ock87gpzRlj87p7hT0EEflaEXNNlvEiyeBvimK+ENmbwAAUcNqpGo01GUZDhn4AJWzP/zX5c27I2veaACgKe/EqChjaHuKCgttNmNoypo2CniMg4NhqwhT2COPwwY++oa7ak6vGL9ObhYJAOh5KVW4j18tlepFrJM1Bhb8z4y2BROtXfbyyQAcu6bDVmcD0t8o60cZH8tw2MBHk2bmBwCeSa2c7qY/JkQq+pmxD2WpewEALpF9LdwTy9M2nEaIrXDIwEdbtqeqWV+EyGXLWUK+poPe4ICyiG4AgKzHe7vx14YQuWhWl2U4ZOCjjbaMDyHGoC/gsdVhFMA+tqOQw5GH/TTRPl6GqQh8DBnqInJQ4KNBNRvEB0HapkgTYgh+VgzsoEjUHAEQ4zijt++ogZux8UXUtljDRhwTBT6Pyc303KEZYEQGe5gVYy8BgymCKkdCAY98NNRlGRT4GEhbwFOasgcAdK4HRGybMael28oQF7+thdyATXWPKDlBhqkDEzlZi5+at6MFDkEZH0PQrC7LoMDnMTE1Ptre8Cjj43gMCXruPi6MNVadSMmJnQAAt+gXjNKeLmK2QxDDWjMrqwr/rPiPhP7Z8iKHDokvY+D/xmt+TRwGBT6EEEKIBTgZONRFIZs8FPgQYkb2ND1cc5jO3K+JipOJLFaU6aGhLsugwMdAmkNcRY8KqOCMmI2ph7h4/OJ9eLyWjbbaJG2L2jED/qhXhQIe68Av12BzwbwVBUDEvCjwkUjXYnT8G4N7ZF+AippJFfhiYTyunbF2/GrFYv5guHt4oLCoyLQdIsQO0Kwuy6DAR6LVBRe0Hh8fNQ4AsJHRcvakai4yA5702SMBAKELPjVmd4zOw91d7WttRd36skGmXhPInM9HrICevRktiQIfy6DAhxAbYu0Bjy7WUNukK9BSDdJsbriG2DSq8bEMhx3cZJyT2k2qO/NGC1PZq5I+e6TwSZ0QRyZ23SLGcWo3U1IdljsW1QXHorqY9PnslbuHh82sSwWgoraH6nsckkNmfIyxH1fN6errsOgaAgOefEofw4Xb5EaVhBgT/+YopQ5IdXVlU87m4hckpMwPMQdnGDjURSO0sjhk4KONOTYnpaCHEN04xvRmdzTvM2UApG2GGtHPLq4VU5q1HsjJwKEuJxNnQ+0V5fke45hSbYNSQoh5yR3SEvM4fedoFmIb2xguXJgNao9s4vXRkBZRQRmfx+RmfPid2/VJbhYJAOh5KVXWcxBijyZ7NAXwZI8nQ2p5DHlsYVGRSYOfNblnTda2NdDMZNt05sfMdT8Gz+qihI8sFPgYmbb9vCjgsX58cBqz+V0A5lsY0JHxb5iWXvNHW9AjZk81sYW8Cu/qsvplqx6smALAePvRmYwVLGBo8KwuGuqShQIfI4n3ag6g4o+5rkUOtbGVdVnsHR+c8pt+EvPRNp1cTAZHs9jZmOdb/Zu2FbOGpQv0omEvh0eBjx5idmzn8bO6ih4V6Ax4ylL3Civg8gJmrjCsk8SoKNNjWfpmfGkWMUsd3tJ1flFhoXBf9vLJAPQHPnxGg3+DJ9q5e3hYb/BjJWioyzIo8NHD0LofzSEv1aCHzwrZ+/g/IXLwWSD+92RV4Z+VzlENhAxd64dvS0ymhwIeYiw01GUZFPiYgGbAo83a4xWZniNPdwdAdUCEaGNNG2Da2v5qRAMNcZHHKPCxEH5IpeclGlohpCr6hk1MvbIzT+7+aoTo4sRxBq3FQ+v4yEOBDyHEJoiZRcXPxqoxbZXOc2hDUvPhv2dlZ/YDeBI8WkMGzxpwzhw4J/nBC0eBjywU+JiArhofbfj9vmrP/cikfSLmkdwskoYtLUiz/ka1SFpuwGPK/aeK8x4CACb5tAFgv6u7W322zEJT252cOTgZEPhQxkceCnxMSExQc/GrUwCAS/9tJ+wTRGwPLVJpn0y96Sa/xo+9Bjy6PFgxhZYMIBZDgY+FlKbsUfuagh7bRgGPbbOpXcXtQMjsDcJwl1Xsi2apwmdnJ3BOBjw3R8O2clCZuwnVnL5SbRd31f1slBFdoYzoii6nDqPLqcPm7xwhhFiQu4eHwwecnBNXUecj9yZjmOzOnTt47bXXULNmTXh6eqJNmzZITX3ywY0xhoSEBISEhMDDwwOxsbG4cOGCWhvFxcWYMGEC/P394eXlhf79++P27dsGXw9zocBHD2NtXFr0qABFjwqwuuACigvyUFyQZ4TeEULEooJm60eBkOllZ2ejU6dOcHV1xb59+3Dx4kWsWLEC1atXF85ZunQpVq5ciQ8//BApKSkICgpCz549kZf35H0rPj4eu3fvxldffYVjx44hPz8f/fr1Q3l5uQVelXQ01KWH3AUMCSGE2BDVD7hG+LArlpMzBycDll92grTHLlmyBKGhodiyZYtwLDw8XPg/YwyrV6/Ge++9h8GDBwMAtm3bhsDAQOzYsQOjR49GTk4OPv30U3z22Wfo0aMHAODzzz9HaGgofvrpJ/Tu3Vv26zEXemfXw1gZHzH44lhCiPFwjInK9pSl7kVZ6l4z9IiQJzgnJ4NvAJCbm6t2Ky4u1vp8e/bsQVRUFF566SUEBASgbdu2+Pjjj4X7r1+/jszMTPTq1Us4plAoEBMTg+PHjwMAUlNTUVpaqnZOSEgIWrRoIZxj7SjwEYEPgIwRBDHOSWcmaQwXrlYHRAghjqKosJDW95EpNDQUvr6+wm3RokVaz/v777+xYcMGNGrUCAcOHMCYMWMwceJEbN++HQCQmZkJAAgMDFR7XGBgoHBfZmYm3NzcUKNGDZ3nWDsa6jKh+8sqNjNUnc7Or/EzqVpLAE82NyWEGB+/qnNVWZ+yiG4A6A+iKr4WUeHlbZbns4rZXWZmrKGu9PR0+Pj4CMcVCoXW85VKJaKiopCYmAgAaNu2LS5cuIANGzbgjTfeEM7TXBiRMVblYolizrEWlPExM3dPL7WFDflM0jOpR7GR3XC49TwIMSWxQ128uwveFlZ/djTJzSKtYsj9WFQXS3fBbAya0fX4BgA+Pj5qN12BT3BwMCIiItSONWvWDLdu3QIABAUFAUClzE1WVpaQBQoKCkJJSQmys7N1nmPtKPCxkDX5aViTnyZ8/Uuk4/yyE0KsT89LqWrrUSm8vI2e7fmpeTv81Lyd3nNoTTPT6dSpEy5fvqx27MqVKwgLCwMA1KtXD0FBQUhOThbuLykpwZEjR9CxY0cAQGRkJFxdXdXOycjIwPnz54VzrB1ldq3EM6lHhfoeyvqQqvBDpjwx26M4IrFDXcS0+J/XzimHRZ2vOa3dXoe/KrI28vMPHKTVnU6ePBkdO3ZEYmIihgwZgt9//x2bNm3Cpk2bKtrjOMTHxyMxMRGNGjVCo0aNkJiYCE9PTwwdOhQA4Ovri5EjR2Lq1KmoWbMm/Pz8MG3aNLRs2VKY5WXtKPAxIdXFC6vCMSUFPES2okcFFPwQq8MHPC6XDlf8G9m3ysfcXfB2pT3X7JW5p7O3a9cOu3fvxqxZszB//nzUq1cPq1evxrBhw4RzZsyYgcLCQowdOxbZ2dlo3749kpKS4O39JPu3atUquLi4YMiQISgsLET37t2xdetWODs7y34t5sQxZt8fhXJzc+Hr64t/Mu6qFX+Zk7Y3JG2LGJqriJDYBn5bE9d2/Svdp5nxASjro43UTMGDFRUTEhzljdfWmCzzozJjNzc3F4HBIcjJyTHZewb/vrQv6ml4ucjPPxSUlaHPqd9N2ld7RBkfM9C2gaXmlHaOKWmoi6gpb94dAOBq4X7YO3cPD+ENlTbOJMT+UeBjJRjnpFbsTIi2DA5fGCq2VoJUcPfweHLtTh21cG8IqeDk7AQnA2p8nBjNT5KDAh8zeCb1yR9afogi3qs5AO3r+Jh7/QxiO/gZL9oCIH1DY45KtUiWv3b2WihLDKCagTfjVkWqU9JlPZ7Zxro51oYCHzOjOgxSFW1Do5r4N3E+2Clv3h3KiK6m75wd0FwojwIhQhwLBT5mNIYLl1S/w2d+7i2fprb6M7FvEa9EqX1dXJBXKft3Z95oAOozB2lTXUJsC2V8LIMCHzNSHdbiAyBts3M0UdDjWDS/39qGPPlzxPz8kIqsDs3YItaGanwsgwIfQohD4GdsiRna4muoaBVhQuwPhYtmVvSoQOundH07v/PDGoRocr7wM5wv/FxRO8Y5AZyTQ+83pUu8Z7MnXzCl2rot2nQ+dZRmfxHTM3SfLgOGyRwZZXwsjN+jq8upwzrP8Z+2XFTBK3E8NINLnI3sxpNMD9VC2ax4z2ZY/eiSpbthNE4cBycnA1ZutpHd0K0NBT4Wws/G4QMZbSs5q+IDI1rk0PEUFRZW2rtI2zmEEEKqRoGPhdCqvMSYVAMjWn1Yu2NRFdlVfUNYYs4hxFg4ZyfDNilVUvZSDgp8CCGEEAsweJNSJQ11yUGBj4Xxhcv+05aLOp+2tXA8VQ1zaeKzFjQjSZ2+lZv5a8xneugaWqfVjy7ZVVbO4HV8KPCRhQIfCxMWoKtilgkhxDg0V24GntTO8YWzFPBYJ5eLB+0i4CGWRYEPIcTh2dNMIXskBKkR3Uz7RCKWOjAmqvGxDAp8bFxpyh6a0kzUULZCHJeLBwEAZRHdhP8ffiMBAF1Da8EHPKrDvSadwfh4LSxzcXKGgTU+RuyMA6HAx8ZR0EOI8VDAY2WoBICYAAU+Nopf90fbPk6EkKq5RPYFAJQVFgr/J5bDr2h/f9kU1JyxGgAqViS3Y5wTB86ABQwNeawjowFCQmxQ+uyRSJ89UtS5fOEu0U7qrDlrNYYLt8nvteY2PjWnr8SxqC7C7C21c1WHucxcj2MKTk5Owkalsm5O9BYuB2V8LEzfHl2E6BK64FNLd4FYGVtbzV3bnoU8vUOO9DeTGIjCRRvHb31BCCG2Qt/Gy/eXTdF5n1p2zsyFyKZgyAalhq4B5Mgo42Oj+NqeoubdadsLopW+T9SEWCthbbOq8JkfGw5+DJ7ObsBjHRkFPjZmUrWWAJ6ktd09vSoVOvNZIJrxZf/EbFpra0MgxP7pC24kFzTbQQBEzIsCHwtjj39ZDan1EbI/jz/hu1PA4zBoCxNiS+RmIUU9jiltLvjhnJzAGVCgbMhjHZlFr1pCQgI4jlO7BQUFCfczxpCQkICQkBB4eHggNjYWFy5csGCPLW9NfprONzt3Ty+7n/5JxLu/bIreeglCrInLpcNwuXTY0t0wK4NmdD2+EeksnvFp3rw5fvrpJ+FrZ2dn4f9Lly7FypUrsXXrVjRu3BgLFixAz549cfnyZXh70/o1hDhdPFTxHy1ZvtpzPzJzbwjRziz1ZpxT5Rlf1p4BMrDGBxT4yGLxwMfFxUUty8NjjGH16tV47733MHjwYADAtm3bEBgYiB07dmD0aN2zAhwBLWBICDEGfoYVHyiLqRszBW2LSJal7q34T7PYqhugae5EJIsHPlevXkVISAgUCgXat2+PxMRE1K9fH9evX0dmZiZ69eolnKtQKBATE4Pjx4/rDHyKi4tRXFwsfJ2bm2vy10CIpVABO5GCDyTKVAIJvtCYz8qsLrig9nVV+OF1U2R1yvhNSaUGNdae6XmMczJwVped1vhMmSJ9iH727Nnw8/MTda5FA5/27dtj+/btaNy4Mf755x8sWLAAHTt2xIULF5CZmQkACAwMVHtMYGAgbt68qbPNRYsWYd68eSbtt6nIKXQuLsiTlPVJbhYJAOh5KVVa5wghNq9MTOZEIlMOY8V7NgOgnn2yp2UaqLhZu9WrVyM6Ohpubm6izj927BjGjx9vG4FPnz59hP+3bNkS0dHRaNCgAbZt24YOHToAADhOfYEmxlilY6pmzZqlFi3m5uYiNDTUyD03PmbAJxTNVLU+XU4dlv08hBiqqLDQbraIsCVSggVtEyQsFWzQUgyOa/fu3QgICBB1rtSaX4sPdany8vJCy5YtcfXqVQwcOBAAkJmZieDgYOGcrKysSlkgVQqFAgqFwtRdtSr+05aLPvfe8mkAqPCVmNfdBW8DAPymilycjliMthofUw5nSWHp5ze2igUMnas+Uefjy43YG+uxZcsW+Pr6ij7/o48+0hsXaLKqPFlxcTEuXbqE4OBg1KtXD0FBQUhOThbuLykpwZEjR9CxY0cL9tL0GOdkUAaIEGsTMnsDQmZvoGyPBfzUvJ2k81cXXMDqgguVNg8FniyZYcyp52KW4SguyENxQR7ivZrr7IstLuXBr9xsyM0eDR8+XFICY+jQofDyEv/9t+hVmzZtGo4cOYLr16/jt99+w4svvojc3FwMHz4cHMchPj4eiYmJ2L17N86fP4+4uDh4enpi6NChluy2yXFMKWtBQ/6Pgz61536E2nM/qvI8QgjRxSWyr9ZZWKai8PKGwssbG9kNceuV2cHO7Y4uOzsb69at0zpBKScnR+d9Ylh0qOv27dt49dVXce/ePdSqVQsdOnTAyZMnERYWBgCYMWMGCgsLMXbsWGRnZ6N9+/ZISkqy2zV8zLlTu8LLm6bEEwCWm75MbJ+7p5fFh59UAzA+I8TPTLN2Tk5OcDKgQNmQx1q7Dz/8EOfOncOECRMq3efr64tffvkFubm5eO+99yS3zTHGmDE6aa1yc3Ph6+uLfzLuwsfHx9Ld0cvYgY/YgIYCIMeVPnskQhd8auluEBPQF5BILWDWl2ExNPAxxRCV0CcZJQO5ubkIDApCTk6Oyd4z+PelC9Neg7dC3MwlbfKKS9B8+ecm7aultGnTBitWrED37t213v/zzz9j2rRpOHPmjOS27TdcJKKHsyZVaylsfkocy+Xd5yzdBWJk2mpzjrWLxbF2sXqHifTdx7f5U/N2lWqG4r2aq9XeiCX3cZJY+ZCXpWt8Fi1aJJSV8MRsFVVcXIwJEybA398fXl5e6N+/P27fvm1QXzT99ddfaNSokc77GzVqhL/++ktW21Y1q4tYBj/EIWVaPLFt/HpOz5w+ZuGeEGPRlXm5v2wKelxIEd2OvhlcnVMOV7pP7rCSKYZW+WFbWxnqsqSUlBRs2rQJrVq1UjsuZquo+Ph4/O9//8NXX32FmjVrYurUqejXrx9SU1PVtp0yhLOzM+7evYu6detqvf/u3buyh/oo8LFzpSl7AABZP1as2KovqLn41anH/0bSAod2jr6/9o+fdeVIH2QqBTxWPjuW4wxcwFDm68vPz8ewYcPw8ccfY8GCBcJxMVtF5eTk4NNPP8Vnn32GHj16AAA+//xzhIaG4qeffkLv3r1lvx5Vbdu2xXfffSes6adp9+7daNu2ray2rfunghhMGdEVyoiuos7tcuowupw6jIhXokzcK0IqPp3zn9CJYbRlZ8w980oKW5x6bgrGGurKzc1Vu6lu26TNuHHj0LdvXyFw4VW1VRQApKamorS0VO2ckJAQtGjRQjjHGMaPH48VK1bgww8/RHn5k/WKysvLsW7dOqxatQrjxo2T1TZlfKyInC0rxBKzyCEtbkjMoaiwEACw+tElta8B7VsUOBKpW9BoY81Bxf1lFavqc85OCJm9wejt86+dH7avOWO10Z/DGmnuTjB37lwkJCRoPferr77C6dOnkZJSefhTzFZRmZmZcHNzQ40aNSqdwz/eGF544QXMmDEDEydOxHvvvYf69euD4zj89ddfyM/Px/Tp0/Hiiy/KapsCHytijuns+mZwUcBDLI0Phog4fKFx5xTbWMDPXH9jLv33NACg8wyzPJ1shhYo849NT09Xm9Wla/G/9PR0TJo0CUlJSXB3d9fdrsStosSeI9XChQsxYMAAfPHFF7h27RoYY+jSpQuGDh2Kp59+Wna7FPg4mKNRsQCA2O1zK+3szafL7y+bQkEQMQr+jZkvrlXN7pDKpGZ7Op86CgA2syK26pCcKQM1W/l5c3J2gpMBgQ//WB8fH1HT2VNTU5GVlYXIyEjhWHl5OY4ePYoPP/wQly9fBqB/q6igoCCUlJQgOztbLeuTlZVlkl0Vnn76aYOCHG2oxsdBubbrL0xR1Zz+yn9aIsRQPS6koMeFFBQVFlr9mxCpoG06vLGZOjtljtdgi7p37460tDScPXtWuEVFRWHYsGE4e/Ys6tevX+VWUZGRkXB1dVU7JyMjA+fPn7eZ7aQo42NFTFnjo42uPz5Spr4S21ac9xAK7+oma58Pdo5FdQEAxG5PQFlEN5M9n6MwxWrb5ggU+Bofc2WU7y+NB2C9tT6cE2fYrC4naUNL3t7eaNGihdoxLy8v1KxZUzjObxXVqFEjNGrUCImJiWpbRfn6+mLkyJGYOnUqatasCT8/P0ybNg0tW7asVCxtrSjwsSLmCHi6nDps8ucgtsOUQQ8APFhR8UYndeiBP89WhnDM7cPTGw16PB/kOF/4GUfiPgDwZI0e5ws/AwDKm2tfMdcQpg54NAPCmtNXmvT5DGWsGh9jErNV1KpVq+Di4oIhQ4agsLAQ3bt3x9atW422ho+p0ZYVVsgcARBtT0FMSTNwMcYsGwqCnig7sx8A4NL2WYPaKXpUIGRh+CDBlIGPyYe4NH7u5AyvmnPLir8S34a3u/hdyDXlFRWjwbsb7HLLClOijI8VMseQ1xgu3GGnDBPTUn2z4f8vfPJmSlGLyvFvXMJKvDTbS42hAY8qzSxMkQkCHnNxPv+47uTxxA0hAHpUYJWLGVpjxscRiAp8+BUcpdi4cSMCAgIkP46Yd5d2QoxF66drmT/Lmmv9kArFeQ8BGG+I0hgZGH1bXIjBP84YfdGcqWrtLLVysz149913kZmZic2bN0t+rKjA57vvvsOQIUPgITLVvGPHDuTn51PgYyBzFzsTIsfdBW8DAPym6q6nUP0ZVvBvlDTLSzJT1mQJw5ES62Kk7I+lLzjSNtVdW1DE/7xpLoBYVFioezjUSgMEztkZTgbUxXA2UlNjCnfu3EF6erqsx4oe6lq7dq3oQObbb7+V1RlSgZnhl3RNfprJn4M4Bn0BjyZDass01wRyNKUpe6wuo+Hu6WWSDUE1AyT+63iv5ngxohYAIGS2xoP0fUBUvc9KgyAizbZt22Q/VlTgc+jQIfj5+YludN++fahdu7bsTpEKlOkh1kxYNfjxInr6qAbzcjM9Yp7HnimbdjFp+9Y8A4rP+Kx+dElYGkEzOJI9VCay7swUqMbHMkQFPjExMZIa7dy5s6zOkAoU8BBbwE9/Fj5Na3vzeHzMGDOy+Dc8R834TPJpY9IJCZqzu8Qw9do/8V7NATyp9XL38HiyNIKI59Y6fKft76vqMTP+/aXAp2pr167VepzjOLi7u6Nhw4bo0qWLpKn0smZ1KZVKXLt2DVlZWVAq1X9IunQx7acSQojljeHChTcjfpE4fnaQakaHX8fHFBtSOhpTz8IUvn8adTVygxt3Ty+jBUa0lIHjWrVqFf799188evQINWrUAGMMDx8+hKenJ6pVq4asrCzUr18fhw4dqrRZqy6SA5+TJ09i6NChuHnzJjSXAOI4Tm37eGIYKm4m1mr1o0vCm5G2gIe/T1vAI3eNFUcf6jI1PjviP225Qe0caxcLwHSZOdWNWW0d52TgrC4DHmsrEhMTsWnTJnzyySdo0KABAODatWsYPXo0/vOf/6BTp0545ZVXMHnyZNH1xZIDnzFjxiAqKgp79+5FcHCw0XdjJYQQQhwBDXVVbfbs2di5c6cQ9ABAw4YNsXz5crzwwgv4+++/sXTpUrzwwgui25Qc+Fy9ehXffvstGjZsKPWhRCLK9BBrI2YrCbHDEvGezQCIW6vH3cODpr+bGJ+5M3QfMFNkelSHzHRleooeFVQqcLbmgm0iTkZGBsrKyiodLysrQ2ZmJgAgJCQEeXl5otuUHPi0b98e165do8DHhGiIi1grYwQ8PP6NVTWg4YMhWlXcwT3+G2iKqfLWhHPiDMv4SNyk1BZ17doVo0ePxieffIK2bdsCAM6cOYO3334b3bpVbHiclpaGevXqiW5TVOBz7tw54f8TJkzA1KlTkZmZiZYtW8LV1VXt3FatWol+ckIIUQ2Y9AU8VOBqHpoZOG1FyqqZFf4+o+68LnfFbzPsMG9MVONTtU8//RSvv/46IiMjhXijrKwM3bt3x6effgoAqFatGlasWCG6TVGBT5s2bcBxnFox84gRI4T/8/dRcTMhhNg2l4sHK/4T2Vc4pmuNHH2BRlnqXqBZrDG7RhxQUFAQkpOT8eeff+LKlStgjKFp06Zo0qSJcE7Xrl0ltSkq8Ll+/bq0nhJCbFbJiZ0AALdo8cWCxH64PA54xNRzqdLM9LhE9hXeYGwtE2MunJMzOCcDtqww4LG2pn79+uA4Dg0aNICLi2H7q4t6dFhYmPD/o0ePomPHjpWeuKysDMePH1c7l1i34oKKYjBDthEg9sfWAh7VNYV4VCtkONWdzTWnqJel7q04qVksXC4drvi/SoaIRwXpVXByrrgZ8ng79+jRI0yYMEHYouLKlSuoX78+Jk6ciJCQELzzzjuS25Q8QNi1a1c8ePCg0vGcnBzJ6SaiHceUZilsVnh5U9BDrNoYLhxjuHBhjRltNrIbcPfwULsR4+qccljWujn0/aiCk5PhNzs3a9Ys/PHHHzh8+DDc3d2F4z169MDXX38tq03J+SK+lkfT/fv34eUlc68UQghRwQ+NaJvVo5lFcLl4UBie4WlmesZw4UJbv0RWrC7f81Kqsbpr11SLm4Uhq8e1O+6eXqg80ZgQ4/nuu+/w9ddfo0OHDmqxR0REBP766y9ZbYoOfAYPHgygopA5Li4OCoVCuK+8vBznzp1Dx44dZXWCEEKAih3IAQDNuwMAnC/8DAAof/w1UHnl57KIbiiTUI/S5dRhY3XXYWgWN4ut2aHaHv04Z2dwEvaY0vZ4e/fvv/8iICCg0vGCggLZCyiLDnx8fX0BVGR8vL294aHyB8bNzQ0dOnTAqFGjZHWCWAbV+BBrkNwsUvj/M6nq21K4tusPAMiaN1r/NOnHQ8P8zzSP/9m29/VgzE01ENLMtpkcvxmuPaxzRjU+VWrXrh327t2LCRMmAIAQ7Hz88ceIjo6W1abowGfLli1gjIExhnXr1sHbm94sbR0FPMSadDl1GEzjmL6MgVrmR8eboRAIcU6V6uYMXaGYaGfyLI89BDxEtEWLFuHZZ5/FxYsXUVZWhjVr1uDChQs4ceIEjhw5IqtNSZVRjDHs2LFDWCaaEEIM1eXUYXQ5dRj3lk/TeY6+bA8Vzxqm6FEBDUlZipPTk6yPrJv9Fzd37NgRv/76Kx49eoQGDRogKSkJgYGBOHHiBCIjI6tuQAtJxc1OTk5o1KgR7t+/j0aNGsl6QqIbbVFBHBEf8JhyXyV3Ty9hSI0fTluTn2ay57MVqgEPv+u5qXZVJ5XRys3itGzZUpjObgySZ3UtXboU06dPx4YNG9CiRQujdYQQ4lg0h5q0ZR10rRis6ZenOgPQX7is+SbB74nniLRda0MDHsoaEVPIyclBcnIybty4AY7jUL9+fXTv3h0+Pj6y25Qc+Lz22mt49OgRWrduDTc3N7UiZwBa1/ghhBApxAY8VeGDm9KUPcIbuyO+QeutlTLwWjvi9TQazsDiZs6+i5s///xzjB8/Hrm5uWrHfX19sXHjRrz88suy2pUc+KxevVrWE5Gq0a7sxJHwmR5+iEXOAnk8YTaYjt8dZQQtrqqNy6XDWldcFotWZjYQzerS6fTp03jzzTcxbNgwTJ48GU2bNgVjDBcvXsTq1avx+uuvo2nTpmjdurXktjmmuvOoHcrNzYWvry/+ybhrUGrMHMwd8NCsLmINDBni4vH1O5pDXaofJjR/3vnnNVZ2yRbo22VdUjt8wGOHH9Jyc3MRGByCnJwck71nCO9Lu9bBx0t+cX5uQSECB08waV8t5c0330R+fj6++eYbrfe/+OKL8PHxwebNmyW3LWunr/Lycnz33Xe4dOkSOI5DREQE+vfvD2cHWEzJVDimpIwPIUam73dJ3x5T9uSn5u3QbMhTAJ4UkEsNeISs3Cn9mTUiDRU36/brr79i/fr1Ou8fM2YMxo4dK6ttyYHPtWvX8Nxzz+HOnTto0qQJGGO4cuUKQkNDsXfvXjRo0EBWRwgFPITIJWc1ZrMvvGdmfPG4MRZvFIYh6W+UcdFQl053795F48aNdd7fuHFj3LlzR1bbkgOfiRMnokGDBjh58iT8/PwAVOzT9dprr2HixInYu3evrI4QQoipOOKwrqEBDxUtmwG/jo8hj7dTjx49UtuUVJNCoUBRUZGstiUHPkeOHFELegCgZs2aWLx4MTp16iSrE4QQAgDxXs0B0GrKhtAXsND1JbbkwIEDwnZZmh4+fCi7XcmBj0KhQF5eXqXj+fn5cHNzk90RQgiR84Y8hgunxQghLkMjJgtUlvo4a/94B3ZiOrRJqX7Dhw/Xe7/cTUol58n69euH//znP/jtt9+EvbtOnjyJMWPGoH///rI6QRx7MTXimIy1VcLqggtQeHmr3czh9txRKC7Iq7QxqiUYc1iqrFksyijoMQ8nJ8NvdkqpVFZ5Ky8vl9W25Ku2du1aNGjQANHR0XB3d4e7uzs6deqEhg0bYs2aNbI64cg4phRumhjnRAERIVUw5aSA0pQ9QnBTnJ+D4vwc4b5aM0y3xYZYfN9oYgQh4kke6qpevTq+//57XL16FX/++ScYY4iIiEDDhg1N0T+7RX+oCFEn93fC6fzPQPuBRu2LkMVRXfjw8YcQa8jwmKIPli5m5r//DvVhj2Z1abVnzx706dMHrq6uos7/8ccf0bVr10o7Segiax0fAGjUqBFtVGpiFBwRR2Doz7mrtqCHb1Pim6g1BDX6iOmfmABCdfFGY6ycbUyqPw/2HgRxTs7gDAheDHmsNRs0aBAyMzNRq1YtUee/8sorOHv2LOrXry/qfMmBT3l5ObZu3Yqff/4ZWVlZUCrV/2gdPHhQapM2TeqnFApmCLE+/Jo3UoukrXmavNi/TdYS8BDTW7RoEXbt2oU///wTHh4e6NixI5YsWYImTZoI5zDGMG/ePGzatAnZ2dlo3749/u///g/NmzcXzikuLsa0adPw5ZdforCwEN27d8f69etRp04do/STMYa4uDgoFApR50ud1i45nJ40aRImTZqE8vJytGjRAq1bt1a72QpddTWmfry11O1MqtbS0l0gpJIxXLgQhBjC6dIROF06Ivr8NflpakHP0ahYHI2K1VswfTQq1uB+mgL/N0bM3xpjXGtTMvTvtNXjDCxslvhecuTIEYwbNw4nT55EcnIyysrK0KtXLxQUPBnmXLp0KVauXIkPP/wQKSkpCAoKQs+ePdVmc8fHx2P37t346quvcOzYMeTn56Nfv36yi401DR8+HAEBAfD19RV1GzZsmKQtOyTv1eXv74/t27fjueeek/xiLEHbXl26ConFkBvsGPJ4U5hUrSWt5UEsih9uUf2d4ANyQ382S1P2AABc24mbacoHANb8OyF2GM4aPljJoe9vo7bXpHm+sV63Offqun/4v/Cp5im/nfxHqBk7RHZf//33XwQEBODIkSPo0qULGGMICQlBfHw8Zs6cCaAiuxMYGIglS5Zg9OjRyMnJQa1atfDZZ58Ju6PfvXsXoaGh+PHHH9G7d2/Zr8dcJP+kuLm52WQhs77ZU6r3G0rqL5+1ZIAIsQaamRdz2chuWHXQI4W9Z0ns/fXJkZubq3YrLi4W9bicnIpZivyCxNevX0dmZiZ69eolnKNQKBATE4Pjx48DAFJTU1FaWqp2TkhICFq0aCGcY+0kv+NOnToVa9asga1t6i4lo6N5U71PzOOl3GepX2J7+SNPiDau7fqLzvYQIuCkDx8ZhN+yQvatoq+hoaFqQz+LFi2q8qkZY5gyZQo6d+6MFi1aAAAyMzMBAIGBgWrnBgYGCvdlZmbCzc0NNWrU0HmOtZNc3Hzs2DEcOnQI+/btQ/PmzStNN9u1a5fROmdMhtbzEEKMi36vxLkzbzQAwH/ackmPs5Xp4YZ+oBT7OKu8DoYuQvj4senp6WpDXWKKgsePH49z587h2LFjle7TXBGZMVblKslizrEWstbxGTRokCn6YpM0f5l0ZXUIIdZJdf0a/nfVmmZr1Z77EQDrn2pv0/i/40xp1h3ojbVlhY+Pj6QanwkTJmDPnj04evSo2kysoKAgABVZneDgYOF4VlaWkAUKCgpCSUkJsrOz1bI+WVlZ6Nixo+zXYk6SA58tW7aIOu/XX39FVFSU6OloxLzGcOE03EUcSnKzSABAz0upasdVP5j8u3QKAKDOvI/N1zEj4T+EmTPTY00ZFY4pK10Doo4xhgkTJmD37t04fPgw6tWrp3Z/vXr1EBQUhOTkZLRt2xYAUFJSgiNHjmDJkiUAgMjISLi6uiI5ORlDhgwBAGRkZOD8+fNYunSpSftfVFSkd8d2sWQvYFiVPn36SFpQyNLoF4YQ23d77igAlQOXMVw41uT9AUB/5sQaAx6xQ126Ah7VgMBYDP07aaq/szb399vMKzePGzcOO3bswPfffw9vb2+hJsfX1xceHh7gOA7x8fFITEwUFilOTEyEp6cnhg4dKpw7cuRITJ06FTVr1oSfnx+mTZuGli1bokePHvJfiw5KpRILFy7Exo0b8c8//+DKlSuoX78+3n//fYSHh2PkyJGS2zRZ4GONxc+mDG5s7heOEDvEPa550JyeviY/Te8QhjUNbWkyRm2P5jF952gyZtBktX8nLdUvMwc+GzZsAADExsaqHd+yZQvi4uIAADNmzEBhYSHGjh0rLGCYlJQEb+8nvyOrVq2Ci4sLhgwZIixguHXrVjibYLf4BQsWYNu2bVi6dClGjRolHG/ZsiVWrVolK/CRvI6PWN7e3vjjjz8snvHh10vIunsb3r7VAVjxL58Z0To+xNK0ZV6MtY6PmOdSZc2BD8+QGh99wYuYv4diPjQa+hymJjaAM+c6Pg9S9sGnmpf8dvIL4Neuj0n7amkNGzbERx99hO7du6vFFX/++Seio6ORnZ0tuU2TZXwIIUQqswQ8Gvt42ULQYwhjZGxsOeCxZpyTk5CllPt4e3fnzh2tawcqlUqUlpbKatPhAx9Hre1ZXXDB0l0gDo7P7phiwcLS376r+E+L7jrPscWAh++zmOyPuYqNJRVTy9w81hjEDPeZHWfgUBdnn5uUqmrevDl++eUXhIWFqR3/5ptvhAJsqUwW+FjjfH4xn1ocLQAixB5oBgL63soU1XxN2xkTcbp4CLDiRRl1BRZVsmAwBFTuJ70HWJe5c+fi9ddfx507d6BUKrFr1y5cvnwZ27dvxw8//CCrTYcqbtbH0X7YHe31EuvDD2vxQYsxMzDK5l0rHbPFDI+qqlaitvR0cp7Wvy1W9PfGqv72cZxhAZ8VJhiM7fnnn8fXX3+NxMREcByHOXPm4KmnnsL//vc/9OzZU1abkgOfwsJCMMbg6VmxsdrNmzexe/duREREqO3dobqTq7VhnJPww29VvwRmZC1/JAkxRkDCt8FPZ681/fFMqMe/304XDgHtBxr8PNbM6oZy9P1tZconb/iO/LfI0C0y7PzalZWVYeHChRgxYgSOHDlitHYlX7UBAwZg+/btAICHDx+iffv2WLFiBQYMGCBMlSOEEEIIMYSLiwuWLVuG8vJyo7YrOfA5ffo0nnnmGQDAt99+i8DAQNy8eRPbt2/H2rVrjdo5Y6Od0O0TvyIvIXXmfay2CKGimi8U1XzhasPZnuKCvEo3qyZ22wczbw9hjfj3JENu9q5Hjx44fPiwUduUPNT16NEjYSGjpKQkDB48GE5OTujQoQNu3rxp1M4RIkbEK1GW7gKxMrZawKxKs/apNGWPUOejGfxoewMU86ZoFUP9UoqbVYfI7AENdVWpT58+mDVrFs6fP4/IyEh4eamve9S/v/SCf8mBT8OGDfHdd99h0KBBOHDgACZPngygYoMya15AyREiY0KIbdOWzRGORVQu2JZLVMBj4dlWWllTX4yB4wwrUHaA4ua3334bALBy5cpK93EcJ2sYTPJP0Zw5czBt2jSEh4fj6aefRnR0NICK7I/cOfWEGKL23I9QnPcQxXkPLd0VQiThh67GcOHCNhtizjcLudkIQ7MYumaFWUN2ipidUqnUeZNb+yM54/Piiy+ic+fOyMjIQOvWrYXj3bt3x6BBg2R1ghBDced/rvhP9AuW7QghImhupqq6YrWUxQlVZ6eadGhL1+NMlYGxxmyTKTg5VdwMeTyRTNY6PkFBQcjPz0dycjK6dOkCDw8PtGvXzioXLSSEEGujaxd4qdkc1QBIzPR1sy7UyjlJy9KYOaNjDYvWGlqg7AglHPPnz9d7/5w5cyS3KTnwuX//PoYMGYJDhw6B4zhcvXoV9evXx1tvvYXq1atjxYoVkjtBiKHcKNNDbFhxfk7Ff0S+kcl9s7aKYmYxzPCGbjPXwsHt3r1b7evS0lJcv34dLi4uaNCggXkCn8mTJ8PV1RW3bt1Cs2bNhOMvv/wyJk+ebLWBj+oPOf3AE0IsSQh0zEDt753mEJLcWVKqjzHHsJS5n89caFZXlc6cOVPpWG5uLuLi4mSX10i+aklJSViyZAnq1KmjdrxRo0YGTWdftGgROI5DfHy8cIwxhoSEBISEhMDDwwOxsbG4cEH+5pqq6WBHF+/V3NJdIMQhmSro0Rw24f/eaR0OUS0WNmXhsKFv7Jo0gzhb/3vOXx9Dbg7Ix8cH8+fPx/vvvy/r8ZKvWkFBgbBdhap79+5BoVDI6kRKSgo2bdqEVq1aqR1funQpVq5ciQ8//BApKSkICgpCz549rXo7DEIIMRdtNSKaH/A4ZZlhQUJVjzPFG7CU/vLn0mwwh/Lw4UPk5Mj7ECF5qKtLly7Yvn07PvjgAwAV8+iVSiWWLVuGrl2lrzORn5+PYcOG4eOPP8aCBQuE44wxrF69Gu+99x4GDx4MANi2bRsCAwOxY8cOjB49WvJzEUKIVZIZOOgdwrenN3x7ei2qaKirSpo7QjDGkJGRgc8++wzPPvusrDYlBz7Lli1DbGwsTp06hZKSEsyYMQMXLlzAgwcP8Ouvv0ruwLhx49C3b1/06NFDLfC5fv06MjMz1TY+VSgUiImJwfHjx3UGPsXFxSguLha+zs3NldwnR7C6QP6QISHEwswRCFhDsGFolsrKAwPGcQbO6rL/mdSrVq1S+9rJyQm1atXC8OHDMWvWLFltSg58IiIicO7cOaxfvx7Ozs4oKCjA4MGDMW7cOAQHB0tq66uvvsLp06eRkpJS6b7MzEwAQGBgoNpxfm8wXRYtWoR58+ZJ6gchhFgtG3kTJ8QUrl+/bvQ2Za/jU9Xc+qqkp6dj0qRJSEpKgru7u87zNNcGYozpXS9o1qxZmDJlivB1bm4uQkNDDeqrvSp6VAAAcPf0quJMQoihtBY16wtmdNWsmIq+AEt1FpicNs3NVoJEGuqq0ogRI7BmzRphj1BeQUEBJkyYgM2bN0tuU9ZV++WXX/Daa6+hY8eOuHPnDgDgs88+w7Fjx0S3kZqaiqysLERGRsLFxQUuLi44cuQI1q5dCxcXFyHTw2d+eFlZWZWyQKoUCgV8fHzUbqQymt1GiI0x1++smIJgay0attZ+6cLv1WXIzc5t27YNhYWFlY4XFhZi+/btstqUHPjs3LkTvXv3hoeHB06fPi3U0+Tl5SExMVF0O927d0daWhrOnj0r3KKiojBs2DCcPXsW9evXR1BQEJKTk4XHlJSU4MiRI+jYsaPUbhNCiPWwtanIUoIhfbOsiDqazq5Tbm4ucnJywBhDXl4ecnNzhVt2djZ+/PFHBAQEyGpb8lDXggULsHHjRrzxxhv46quvhOMdO3aUNPzl7e2NFi1aqB3z8vJCzZo1hePx8fFITExEo0aN0KhRIyQmJsLT0xNDhw6V2m2Ts4blzwkhNsCe3qys6e+d1C0yiFWrXr06OI4Dx3Fo3Lhxpfs5jpNdzys58Ll8+TK6dOlS6biPjw8ePnwoqxO6zJgxA4WFhRg7diyys7PRvn17JCUlVRrrswYU8BBC9LKngMfa2UhBOO3VpduhQ4fAGEO3bt2wc+dO+Pn5Cfe5ubkhLCwMISEhstqWHPgEBwfj2rVrCA8PVzt+7Ngx1K9fX1YneIcPH1b7muM4JCQkICEhwaB2SWX2/AtDiDWRuvEokUlfMbi1/r3jDNyd3VpflxHExMQAqJjVFRoaCicj7kQvOfAZPXo0Jk2ahM2bN4PjONy9excnTpzAtGnTZG0WRmxfcd5DAIDCu7pF+0GItSh6VCBvtqTq3llyZ1KRyrTtSaZ6XVWvNV1vqxIWFgYAePToEW7duoWSkhK1+zV3fBBDcuAzY8YM5OTkoGvXrigqKkKXLl2gUCgwbdo0jB8/XnIHCCHErll71sHRmHuZAH1oOnuV/v33X7z55pvYt2+f1vvLy8sltynpqpWXl+PIkSOYOnUq7t27h99//x0nT57Ev//+K2xhQRyPwrs6ZXsIUfFLZBcUF+RVDHNJeXNzgDcyi9HM5ljDrCia1VWl+Ph4ZGdn4+TJk/Dw8MD+/fuxbds2NGrUCHv27JHVpqSMj7OzM3r37o1Lly7Bz88PUVFRsp6UWJfSlIofHtd2/S3cE0LsQ89LqU9qe8RkfLSdQ0MupkXX1yYcPHgQ33//Pdq1awcnJyeEhYWhZ8+e8PHxwaJFi9C3b1/JbUoOF1u2bIm///5b8hMR61XevDvKm3e3dDcIsS98hkHfJ3OqKdHP3rMalPGpUkFBgbBej5+fH/79918AFbHI6dOnZbUp+aotXLgQ06ZNww8//ICMjAy1RYVoQ1DbwjElTcMnxBJ0BTz0++hQ+E1K5d/sf+XmJk2a4PLlywCANm3a4KOPPsKdO3ewceNGyfuD8iQXN/PbwPfv319tzyx+Dy05hUbEsvjgh/buIsTCHOATvCQUCDq8+Ph4ZGRkAADmzp2L3r1744svvoCbmxu2bt0qq03Jgc+hQ4dkPREhhBBiEKmbulo7mtVVpWHDhgn/b9u2LW7cuIE///wTdevWhb+/v6w2JQc+/KJCxLZobqmhOsRFixkSYhh+53VFNd8nB/X9Xulao8cW37yJfIZuNGrnQ12lpaVo0qQJfvjhB0RERAAAPD098dRTTxnUruTA59y5c1qPcxwHd3d31K1bFwqFwqBOEePTVstDAQ8hppHcLBJdTh22dDeItaOMj16urq4oLi5WK6sxBslXrU2bNmjbtm2lW5s2bdC0aVP4+vpi+PDhKCoqMmpHiXGoFsbpwtf6EELk6ZJy0NJdsE+UETOK9evXo169enB3d0dkZCR++eUXS3dJpwkTJmDJkiUoKyszWpuSA5/du3ejUaNG2LRpE86ePYszZ85g06ZNaNKkCXbs2IFPP/0UBw8exOzZs43WSUIIIQAYU78Rm2bYjC55G5x+/fXXiI+Px3vvvYczZ87gmWeeQZ8+fXDr1i0TvELD/fbbb9i1axfq1q2L3r17Y/DgwWo3OSQPdS1cuBBr1qxB7969hWOtWrVCnTp18P777+P333+Hl5cXpk6diuXLl8vqFDEdqu0hxHTGcOEAgDV5f1TOTmj7fZOyHxcFOhXsKetjgaGulStXYuTIkXjrrbcAAKtXr8aBAwewYcMGLFq0SH5fTKR69ep44YUXjNqm5MAnLS1N2DRMVVhYGNLS0gBUDIfx08+I5TDOSdQ6Pfw5FAgRIo9aUTPPHL9PqsGQnRe6yuIgK2FrrqGnUCi01tqWlJQgNTUV77zzjtrxXr164fjx4ybto1xbtmwxepuSfzObNm2KxYsXq+2QWlpaisWLF6Np06YAgDt37iAwMNB4vSSyVBX0aC5gqPr1GC5c+PRKCBFnTX4a1uSnmT7o4WcDqd6IzalYwNCwGwCEhobC19dXuOnK3Ny7dw/l5eWV3p8DAwORmZlp8tcrV1lZGX766Sd89NFHyMur2Arm7t27yM/Pl9We5IzP//3f/6F///6oU6cOWrVqBY7jcO7cOZSXl+OHH34AAPz9998YO3asrA4RYmzps0cCAEIXfGrhnhBiIApw7IqhpVr8Y9PT0+Hj4yMcr2pmteYsKX4BYmt08+ZNPPvss7h16xaKi4vRs2dPeHt7Y+nSpSgqKsLGjRsltyk58OnYsSNu3LiBzz//HFeuXAFjDC+++CKGDh0Kb29vAMDrr78uuSPEvKoaBltdcMGMvTEtCniI1bLj4RdiPj4+PmqBjy7+/v5wdnaulN3Jysqy2lGaSZMmISoqCn/88Qdq1qwpHB80aJBQpySV5MAHAKpVq4YxY8bIekJCCLE3knZiVz1PDCpqlk/b90NKQbmJKRmD0oDvr9THurm5ITIyEsnJyRg0aJBwPDk5GQMGDJDdD1M6duwYfv31V7i5uakdDwsLw507d2S1KWsg+rPPPkPnzp0REhKCmzdvAgBWrVqF77//XlYniPmJ3Zy06FEBretDSBUUXt5QeHk7zI7ZNskKvzfMCDeppkyZgk8++QSbN2/GpUuXMHnyZNy6dctqkxlKpVLrHqC3b98WRpmkkvxTsGHDBkyZMgV9+vRBdna20KEaNWpg9erVsjpBCCF2gd913ZjZBCpgNi5jf39szMsvv4zVq1dj/vz5aNOmDY4ePYoff/xR62xta9CzZ0+12ILjOOTn52Pu3Ll47rnnZLXJMSYtVxYREYHExEQMHDgQ3t7e+OOPP1C/fn2cP38esbGxuHfvnqyOmEpubi58fX2Rdfe2qDFQUkFzajvt2E5I1fg9uwBUzi7IfbPl/0RT4COftkyPju9Hbm4uAmrXRU5OjsneM/j3pVt3Mw16jtzcXNQNCTJpXy3t7t276Nq1K5ydnXH16lVERUXh6tWr8Pf3x9GjRxEQECC5Tck1PtevX0fbtm0rHVcoFCgooCERe0Fr+xBiIGNnFRij4MfOMMYgMfdQ6fH2LiQkBGfPnsWXX36J06dPQ6lUYuTIkRg2bBg8PDxktSk58KlXrx7Onj1bKS22b98+YfdUQgghRkZBj2GscHhLySpuhjzeEXh4eGDEiBEYMWKEUdqTHPhMnz4d48aNQ1FRERhj+P333/Hll19i0aJF+OSTT4zSKUIIIYQQALh8+TLWrVuHS5cugeM4NG3aFOPHjxcWTZZKcuDz5ptvoqysDDNmzMCjR48wdOhQ1K5dG2vWrMErr7wiqxPEeomd/UUIMTKq7TEeK/475iBJG9m+/fZbvPrqq4iKikJ0dDQA4OTJk2jZsiV27NiBl156SXKbkoubVd27dw9KpVJWcZG5UHGzcSi85E0bJMSRqBU3E5tkzuLmv9Iz4G3Ac+Tl5qJBaLBdFzfXr18fr732GubPn692fO7cufjss8/w999/S27ToMpVf39/qw56iOEY5wTGOaE0ZQ9KU/ZYujuEWLVJ3q0xybu1/pMcoCCVEGPJzMzEG2+8Uen4a6+9Jnt/MVFDXW3bthW9j8fp06dldcTkrDjVac34oa7y5t0BAK6W7Awhtkw14NE1jEW7rTsUmtVVtdjYWPzyyy9o2LCh2vFjx47hmWeekdWmqMBn4MCBwv+Lioqwfv16REREqI23XbhwgTYmJYQ4tHUnVwMA1D5miXlzooDHISmh8bMi4/H2rn///pg5cyZSU1PRoUMHABUxxzfffIN58+Zhz549aueKIbnG56233kJwcDA++OADteNz585Feno6Nm/eLKU5kxNqfO7cgo9vdUt3x2bx6/nQQoaE6Fb623cAAGXzrk8O6vsTywc5VMhsfjr27DJnjc+VW3cNrvFpXDfErmt8nJzEVeRwHKd1awutbUrtxDfffKNzvG3nzp1Sm7M81SXmaThMJ44pwTEl7dtFiCnQlhQOiTHDb/ZOqVSKuokNegAZgY+HhweOHTtW6fixY8fg7u4utTlCCLF5xXkPUZz3UP2go7wz2Sor+LDLL2BoyI1IJ3kdn/j4eLz99tuVxts2b96MOXPmGL2DJkdbMhBCZOKHtvB4aMu1/UAANK2dEGP6/fffcfjwYWRlZUGpVA9WV65cKbk9yYHPO++8g/r162PNmjXYsWMHAKBZs2bYunUrhgwZIrkDxPYkN4sEAPS8lGrhnhBiWXwtj6KaLwCVgKeqTA8NaxHQrC4xEhMTMXv2bDRp0gSBgYFqM8zFzjbXZNAChragygUM+VQnZX700rZZKRU6E1KhUoZHTEEzsUrmLG4+f/2OwcXNLerVtuvi5sDAQCxZsgRxcXFGa1Nyxsem6QtyKADSi3ZrJ0QPMZ8fKeAhGhgMKwOz66zFY05OTujUqZNx2xRzkp+fH+7duye60bp16+LmzZuyO2UWVlDYRgghhBDdJk+ejP/7v/8zapuiMj4PHz7Evn374OvrK6rR+/fvS5paZlWYkrI+hJAqFefnwOnCoYovImIt2hdim5SMQWlAyseQx9qKadOmoW/fvmjQoAEiIiLg6qq+f8CuXbsktyl6qGv48OGSG7daqsGNZtaHgh69OKak4S7i0FTreYSFCh3gDYgYH4Nhw1WO8FM3YcIEHDp0CF27dkXNmjVlFzSrEhX4aE4fs0m6hrU0AyDK+IhWXJCn9rXTxUNwbSduyXBCbM0YLhwAsCbvD8t2hBAHsn37duzcuRN9+/Y1WpuOVdysiup7jE4Z0bXqkwixUbIDHipqJjoYugihIyxg6OfnhwYNGhi1TUptaEOFz6JMqtYSk6q1VDt2Z95o3Jk32kI9IsTMaHVmYghDt6twgB+9hIQEzJ07F48ePTJam46b8SGycRpBocLLW/i//7Tl5u4OIZZBAQ8hJrd27Vr89ddfCAwMRHh4eKXi5tOnT0tukwIfbajGR5TVBRcAVK71UT3GZ4Q2shtm6xchVoGGuEgVlGBQGpC2MeSxtmLgwIFGb5MCH210LWaomumg4IgQu6d1zy3K9BAjMXSk1BF+FOfOnWv0NmW9e//111+YPXs2Xn31VWRlZQEA9u/fjwsXLhi1c8S6cUxZadhL05r8NKzJTzNTjwgxjuL8HNpolBAr8fDhQ3zyySeYNWsWHjx4AKBiiOvOnTuy2pMc+Bw5cgQtW7bEb7/9hl27diE/Px8AcO7cOZNEZlaBip1lczr/M5zO/0xvIsS2qH4Ur1RRSohx8LO6DLnZu3PnzqFx48ZYsmQJli9fjocPHwIAdu/ejVmzZslqU3Lg884772DBggVITk6Gm5ubcLxr1644ceKErE5YLW0BD+dEw1wSuLYfCNf2Ay3dDUKk4TjDanSovoeIYMiMLkeJw6dMmYK4uDhcvXoV7u7uwvE+ffrg6NGjstqU/A6elpaGQYMGVTpeq1Yt3L9/X1YnrJZmkMM5UfZHD6eLhyodKy7Iqyh05pxo+IBYveRmkUhuFmnpbhAHwRc3G3KzdykpKRg9uvISKbVr10ZmZqasNiUHPtWrV0dGRkal42fOnEHt2rVldcJqSQlyKCCiBQwJ4TnCR3FCzMDd3R25ubmVjl++fBm1atWS1abkwGfo0KGYOXMmMjMzwXEclEolfv31V0ybNg1vvPGGrE7YDNXAhg90KODRi6/xIcQWdEk5iC4pB58coACGmBANdel269YtKJVKDBgwAPPnz0dpaSkAgOM43Lp1C++88w5eeOEFWW1LDnwWLlyIunXronbt2sjPz0dERAS6dOmCjh07Yvbs2bI6YdP44TCq/dFKrcaHrhGxUmO4cGEvLrNwlHcvohe/O7shN3tVr1493Lt3D8uXL8e///6LgIAAFBYWIiYmBg0bNoS3tzcWLlwoq23J6/i4urriiy++wPz583HmzBkolUq0bdsWjRo1ktUBq6Frt3ZiEG2LGxJibSrtw2XHbyiE2AL2+HfQx8cHx44dw8GDB3H69GkolUo89dRT6NGjh+y2ZS9g2KBBA6NvHGZTxGQudC2E6KgoqCSOgmZ1ERHKlRU3Qx7vKLp164Zu3boZpS1Rgc+UKVNEN7hy5UrZnbEoelM2vccB4O25owAAdeZ9bMneEAJAY3VmyvQQMzJ0uMqeh7oA4JNPPkG1atX0njNx4kTJ7YoKfM6cOaP2dWpqKsrLy9GkSRMAwJUrV+Ds7IzISAeaBqptfR9CiM0o/e27iv80N8FsRMZsLutztF3Fp2m14m7ikG7cuIEPPvgABw8eRGZmJkJCQvDaa6/hvffeU1u/79atWxg3bhwOHjwIDw8PDB06FMuXL1c7Jy0tDePHj8fvv/8OPz8/jB49Gu+//z44Eb8fGzduhLOzs877OY4zXeBz6NCT9VlWrlwJb29vbNu2DTVq1AAAZGdn480338QzzzwjuQN2hykpCKpCrRk2mhUkdkWpLeDh/xib4pO0lX8673kpFYCO/cmISSgZQ7kVZnz+/PNPKJVKfPTRR2jYsCHOnz+PUaNGoaCgAMuXLwcAlJeXo2/fvqhVqxaOHTuG+/fvY/jw4WCMYd26dQCA3Nxc9OzZE127dkVKSgquXLmCuLg4eHl5YerUqVX249SpUwgICDD66+MYk3blateujaSkJDRv3lzt+Pnz59GrVy/cvXvXqB00VG5uLnx9fZF15xZ8fHzM86S6CqUpIFKj8PK2dBeIAzPpBqTaPs3qa9vGskP2LDc3FwG16yInJ8dk7xn8+1LSuevw8pb/d7AgLw+9WtUzaV95y5Ytw4YNG/D3338DAPbt24d+/fohPT0dISEhAICvvvoKcXFxyMrKgo+PDzZs2IBZs2bhn3/+gUKhAAAsXrwY69atw+3bt/VmfZydnZGRkWGSwEfyO3Fubi7++eefSsezsrKQl0czePQGN7TuDyEWd3vuqIo6M5pOTohoOTk58PPzE74+ceIEWrRoIQQ9ANC7d28UFxcjNTVVOCcmJkYIevhz7t69ixs3buh9Pok5GUkkz+oaNGgQ3nzzTaxYsQIdOnQAAJw8eRLTp0/H4MGDjd5Bu0NZHwG/NQCfYifEHGpNW2b6J7HBGh9ifsaa1aW5srFCoVALNgz1119/Yd26dVixYoVwLDMzE4GBgWrn1ahRA25ubsJWEpmZmQgPD1c7h39MZmYm6tWrp/M5586dW2Vhs1yS34U3btyIvn374rXXXkNYWBjCwsIwbNgw9OnTB+vXrzdFHwkhxLQsnfmx9PMTizDWAoahoaHw9fUVbosWLdL6fAkJCeA4Tu/t1KlTao+5e/cunn32Wbz00kt466231O7TNlTFGFM7rnkOn8mpqrh57ty58PT01HuOXJIzPp6enli/fj2WLVuGv/76C4wxNGzYEF5eXqbon20SM5RFa/ygy6nDos8dw4VjI7thsr4QYhUoU+RQyg0sbuYfm56erlbjoyvbM378eLzyyit621TN0Ny9exddu3ZFdHQ0Nm3apHZeUFAQfvvtN7Vj2dnZKC0tFbI6QUFBlTYSzcrKAoBK2SJzkr2AoZeXF1q1amXMvtg2sQGMZlBEAZBe6bNHAgDW5J61bEeIzeOHVrv8bua94yiTQ0zMx8dHVHGzv78//P39RbV5584ddO3aFZGRkdiyZQucnNTfo6Kjo7Fw4UJkZGQgODgYAJCUlASFQiEsbRMdHY13330XJSUlwhT3pKQkhISEVBoCMyfJgU/Xrl31pqgOHqQ1IIhpFOc9BAAovKtbtB+EVInj5Ac8/OMo82P3lACUBsTFppomc/fuXcTGxqJu3brCXlm8oKAgAECvXr0QERGB119/HcuWLcODBw8wbdo0jBo1SgjChg4dinnz5iEuLg7vvvsurl69isTERMyZM0fUOj6mIjnwadOmjdrXpaWlOHv2LM6fP4/hw4cbq1+2w9BMDWV6hM0hqxzKevyLwi88J2x+Sog1MnRNIAqA7F65kqHcgMjHkMfqk5SUhGvXruHatWuoU6eO2n18jY6zszP27t2LsWPHolOnTmoLGPJ8fX2RnJyMcePGISoqCjVq1MCUKVMk7QZRWFgIxphQ73Pz5k3s3r0bERER6NWrl6zXJ3kdH10SEhKQn5+v9qKtgcnX8dEMXKROVVdd88fBgyBt6/rw21vUmras8hvJ468V1XzN0j9i2/iModkYazFECnzMypzr+Ow8dQ1e1QxYxyc/Dy9ENTTLOj6W0qtXLwwePBhjxozBw4cP0bRpU7i6uuLevXtYuXIl3n77bcltGu2d9rXXXsPmzZuN1RwhhBBi15iBM7pMudaNtTh9+rSwK8S3336LwMBA3Lx5E9u3b8fatWtltSm7uFnTiRMn4O7ubqzmbIehixHSYoaYVK0lAO1DXfxGpmqr7D7+BPzvsmlq5xCilym3o9DGWM9DQ152q5xV3Ax5vL179OgRvB+vbp2UlITBgwfDyckJHTp0wM2bN2W1KTnw0VykkDGGjIwMnDp1Cu+//76sTpDHHHSG15r8tKpPomm+xFAO8OmYEHvTsGFDfPfddxg0aBAOHDiAyZMnA4CwLYYckgMfHx8ftWpsJycnNGnSBPPnz5ddaEQIIWZj7swPITqoLkIo9/H2bs6cORg6dCgmT56M7t27Izo6GkBF9qdt27ay2pQc+GzdulXWExFibLWmW1chPRGHH7a0WFG6sYqNHeBNh5iWtc7qsiYvvvgiOnfujIyMDLRu3Vo43r17dwwaNEhWm5LHVOrXr4/79+9XOv7w4UPUr19fUlsbNmxAq1athMWXoqOjsW/fPuF+xhgSEhIQEhICDw8PxMbG4sKFC1K7TAghhBAbFRQUhLZt26otovj000+jadOmstqTnPG5ceMGysvLKx0vLi7GnTt3JLVVp04dLF68GA0bNgQAbNu2DQMGDMCZM2fQvHlzLF26FCtXrsTWrVvRuHFjLFiwAD179sTly5eFYifiIKi+h1gLyvQQI6GhLu0GDx6MrVu3wsfHp8rNz3ft2iW5fdGBz549e4T/HzhwAL6+T9LU5eXl+PnnnyUvQf3888+rfb1w4UJs2LABJ0+eREREBFavXo333ntPeOHbtm1DYGAgduzYgdGjR0t6LpuiOdPLQYqdiwvytK7lQ2wXvw7ThS9TAQBdt70HRfQLaudoW1/HZlbnNsewFwX9dotmdWnn6+sr1BKrxhrGIjrwGThwIICKHVU1V2h2dXVFeHi42pb1UpWXl+Obb75BQUEBoqOjcf36dWRmZqoVTCsUCsTExOD48eM6A5/i4mIUFxcLX+fm5sruE7ENyc0i0fNSqqW7QR4rznuIST5tAACDGvkBAGJTK7ayYRC3kCBtT6IDTW23K5Tx0W7Lli1a/28solMJSqUSSqUSdevWRVZWlvC1UqlEcXExLl++jH79+knuQFpaGqpVqwaFQoExY8YIS1HzO7pq7uAaGBhYabdXVYsWLYKvr69wCw0NldwnYn2cLhyydBeISArv6tjIbqity5S1ZCqylkyV3FbJiZ0oObFT633FeQ/NvxqzJsZMP/TFP4edvskRok9hYSEePXokfH3z5k2sXr0aSUlJstuUXONz/fp12U+mTZMmTXD27Fk8fPgQO3fuxPDhw3HkyBHhfs2NzBhjejc3mzVrlto+ILm5ubYf/NB2FnBtP1B9EUNiE/hMj1ysRXcj9QQ4+nRFW2bfnd1Q2obTKPNjF5RKBqUBM7MMeaytGDBggNqWFU8//TTc3NwM2rJCVOCzdu1a/Oc//4G7u3uVS0RPnDhRUgfc3NyE4uaoqCikpKRgzZo1mDlzJgAgMzNT2PIeqFi0SDMLpEqhUEChUEjqA7FtXVIMe3MlpnM4shsAIGJIxXobATPlD4fzLJ7lMTd9mR5t91EwZDOUBtb4OEDcg9OnT2PVqlUAKrasCAoKwpkzZ7Bz507MmTPHdIHPqlWrMGzYMLi7uwsd0IbjOMmBjybGGIqLi1GvXj0EBQUhOTlZWKSopKQER44cwZIlSwx6DpvkIKs6l6ZUFNG7tutv4Z44HlMUGRua8eGlzx4p/F8zeCrOe0i1QKooG0TsiMW2rFAd3jLmUNe7776LPn36IDQ0FHl5efjqq69w+PBh7N+/HxzHIT4+HomJiWjUqBEaNWqExMREeHp6YujQoUbrA7EuyoiuOu/TXPBOdehrDBcOQPt+X0Q/oYZGy7CSruwKd/5nuGnMztKGD0gMzdIYI1PEF8A7XMaIWC0qbq6aKbaskJw+mD9/vlqhEa+wsBDz58+X1NY///yD119/HU2aNEH37t3x22+/Yf/+/ejZsycAYMaMGYiPj8fYsWMRFRWFO3fuICkpyTHX8OGc7D7bo+rOvNG4M0/7zL3S375D6W/fmbdDdoy16C65lkbs+QcaP4UDjZ+S0y1iKCqItnrljBl8s3dz5szBtGnTEB4ejvbt2xtlywqOSdzX3tnZGRkZGQgICFA7fv/+fQQEBGhd3NCScnNz4evri6w7t2RHh1bBgYIeVdrW9eGzOy80rQmgoliVnz5NGR9p0mePlJ1NkTK8ZOosCz9jLHTBpxbvi9WhIS9JcnNzEVC7LnJyckz2nsG/L3146Dw8qsn/IF+Yn4fxXVuYtK/WIDMzU9iygl+9+ffff4ePj4+s1Zslz+rSNavqjz/+gJ+fn+QOEPs2qVpLACJ3YNciuVkkAKit06MZ3BTnPaSARyZDhpA019qx5EKExhgKs1uMWS74oXojvWhWV9W2bt2Kl19+GUFBQWrHn376adltig58atSoAY7jwHEcGjdurBb8lJeXIz8/H2PGjJHdEUK06XLqsM77+FWBa01bZqbe2A9++MkYxcf62uKLks0VmGgLlDX9u3w6APq5MYiYgMYBhmEMVQ4DV242Wk+s16xZszBx4kS89NJLGDlyJDp27Ghwm6KHurZt2wbGGEaMGIHVq1erLSPt5uaG8PBwYezNmtBQl32YVK1l5UwPX9zMGM3qkcjeh3vE/DzY+zXQyliZF32Bj+pbig1mesw51LXq5zR4GLBNT2FBHiZ3b2nXQ13l5eXYu3cvtm7dir1796JevXp48803MXz48EpZILFEZ3z4bSrq1auHjh07wtXVVdYTEiNxkOnt+hxtV7FGjOqCdLTVgX58fdSa3LMW7QcxM80AROoQFGVvTIJmdVXN2dkZ/fv3R//+/ZGVlYXPP/8cW7duxfvvv49nn30WI0eOxPPPP6+2c3tVJNf4xMTECP8vLCxEaWmp2v32GnUSy9JWI6Q6nOGQn95loIDHwcl5oxT7GAd4EzY2Q2dmOcKsLlUBAQHo1KkTLl++jCtXriAtLQ1xcXGoXr06tmzZgtjYWFHtSE4XPHr0COPHj0dAQACqVauGGjVqqN0IMTcKeqpmFftaWRnu/M/gztvY9hVySZnarro3mIO9sZqbUslQbsDNEYqbgYqlb5YvX47mzZsjNjYWubm5+OGHH3D9+nXcvXsXgwcPrrR5uj6SA5/p06fj4MGDWL9+PRQKBT755BPMmzcPISEh2L59u9TmiFhMqX5TPa7vMbra0fy/lTkaFYujUbGiztW2+aW+zS2J/Tr6dHdhT66qyFm/yO4YGthYe4BkzX0jojz//PMIDQ3F1q1bMWrUKNy5cwdffvklevToAQDw8PDA1KlTkZ6eLrpNyUNd//vf/7B9+3bExsZixIgReOaZZ9CwYUOEhYXhiy++wLBhw6Q2SUxFtf5HM8DRFjxZUb2QsMpuQZ6sxzv8G5qKkhM7ta7KTIhajQ8FCGbHZ24Meby9CwgIwJEjR/ROngoODpa0q4TkwOfBgweoV68egIp6ngcPHgAAOnfuLGuzMGICqoGMlWZ0jCl0wac0jKMHBYFEJ1sLdqSuSWTls8oo8Knap59WvSgpx3EICwsT3abkwKd+/fq4ceMGwsLCEBERgf/+9794+umn8b///Q/Vq1eX2hwxBl0ZG7FBjxVlenj8eiz61vERI332SFGr+dojc6+hY2ucLh4GACgjYi3aDyKCrQVoxCBr164Vfa6cjdElBz5vvvkm/vjjD8TExGDWrFno27cv1q1bh7KyMqxcuVJyB4gRyc3uWOFQFx/w/Lt0CgBgwfxkWaszO/KbviO+9mZDaF8wq2Loys3GCngsuXq1HuVKw7I25Xaa0F+1apWo8ziOkxX4SN6rS9OtW7dw6tQpNGjQAK1btzakKZOwmwUMzcGKAh+etr26tBEz1OVoa/s44vDf4ciKtZ16Xzkt+jGOeJ0sRmrwoe3tiW/DRNthmHMBw9nfp8Ldq5rsdooK8rFgQKRdL2BoCga/09WtWxeDBw+Gn58fRowYYYw+EXOzg53fHWpqchWSm0UKQ4WEWBWaZUWsgOShLl0ePHiAbdu2YfPmzcZqkhBRey8BgFv0CwDo07ujc3K2vuEMIpLUDI621aitcDhLHypuFuf27dvYs2cPbt26hZKSErX75JTYGC3wIQ6CKc2aHTK0uNkRqW7h4Wgc+bUT26M0MPBxhAUMf/75Z/Tv3x/16tXD5cuX0aJFC9y4cQOMMTz1lLyaPtse3yCWYcULHxJCbAw//MVxFbeqhsJ0DZfZWLaHiDNr1ixMnToV58+fh7u7O3bu3In09HTExMTgpZdektUmBT5EWiBjoVogsQsZKryr6yxipm0bCLEBfPCjKwDiAySeDdcNlTPDtqxwhL26Ll26JGxH4eLigsLCQlSrVg3z58/HkiVLZLUpeqhr8ODBeu9/+PChrA4QG2MHmR7awd3+8bvQa1sCgb7/Ns6O3uypxqdqXl5eKC4uBgCEhITgr7/+QvPmzQEA9+7dk9Wm6MDH19e3yvvfeOMNWZ0gRAw+6yN2ijshxqAZIFHW0Mg0p6c7EAp8qtahQwf8+uuviIiIQN++fTF16lSkpaVh165d6NChg6w2RQc+W7ZskfUExAbww1fatrrQvE+VFS58yJvk0wYAsCb3rEX7QSxD3/ddM5BReFfXGcxw538GHs8YFDa9pS1ADKc6A8sBAx5bUlxcjPbt2+OPP/7AmTNn0KZNG+G+W7duYdy4cTh48CA8PDwwdOhQLF++HG5ubsI5aWlpGD9+PH7//Xf4+flh9OjReP/998GJrMlauXIl8vPzAQAJCQnIz8/H119/jYYNG4pe6FATzeoilUkdzhITMBkxOOIzP/yqznXmfSyvnbyHdjncwe9OTjOc5BN+Lh4HPcCTPc/o+hqJuQMeEy14aIgyJYOzAVmbMjNkfGbMmIGQkBD88ccfasfLy8vRt29f1KpVC8eOHcP9+/cxfPhwMMawbt06ABULNfbs2RNdu3ZFSkoKrly5gri4OHh5eWHq1Kminr9+/frC/z09PbF+/XqDXxMFPkR/oCMlCNK247ucdgBRgVKtGbrXb+BrO2hYgsgh7BWnJbhp+kIbM/eGiCJmlWcrY+1DXfv27UNSUhJ27tyJffv2qd2XlJSEixcvIj09HSEhIQCAFStWIC4uDgsXLoSPjw+++OILFBUVYevWrVAoFGjRogWuXLmClStXYsqUKaKzPgBQUlKCrKwsKJXq7yV169aV/Loo8CHGY8zd4CVkipKbRVa5wCFxTGIKmbUFx/qyOfymtxRUWwl9wY2+Y1aU+TFUbm6u2tcKhQIKhcKgNv/55x+MGjUK3333HTw9PSvdf+LECbRo0UIIegCgd+/eKC4uRmpqKrp27YoTJ04gJiZGrS+9e/fGrFmzcOPGDdSrV6/Kfly5cgUjR47E8ePH1Y4zxsBxHMrLyyW/Ngp8iHHpGuIy4WwwWuRQnbLcOj/dEkLUGWsBw9DQULXjc+fORUJCgux2GWOIi4vDmDFjEBUVhRs3blQ6JzMzE4GBgWrHatSoATc3N2RmZgrnhIeHq53DPyYzM1NU4PPmm2/CxcUFP/zwA4KDgyVliXShwIcYj74CaBPTN33Z0fAbdFJGwjQONK5YLTY29aCFe0IAiB/G0swMWUHmp5wZthYP/9j09HS1TUp1ZXsSEhIwb948vW2mpKTg+PHjyM3NxaxZs/Seqy0I4TMxus7h90UXG8CcPXsWqampaNq0qajzxaDAh1g3kUNea/LTtB7/d/l01Jq2zNi9IoTYIxtdDNHHx0fU7uzjx4/HK6+8ovec8PBwLFiwACdPnqwUQEVFRWHYsGHYtm0bgoKC8Ntvv6ndn52djdLSUiGrExQUJGR/eFlZWQBQKVukS0REhOz1enShwIfYtQtfpiJ2mqV7QQixOH0BjYWCHXMXN/v7+8Pf37/K89auXYsFCxYIX9+9exe9e/fG119/jfbt2wMAoqOjsXDhQmRkZCA4OBhARcGzQqFAZGSkcM67776LkpISYYp7UlISQkJCKg2B6bJkyRLMmDEDiYmJaNmyJVxdXdXuFxPwaaLAhxBCCLEAa53VpTlTqlq1agCABg0aoE6dOgCAXr16ISIiAq+//jqWLVuGBw8eYNq0aRg1apQQjAwdOhTz5s1DXFwc3n33XVy9ehWJiYmYM2eO6KGuHj16AAC6d1dfP4uKmwkhRA9jrtlEtT12xgaHtqyBs7Mz9u7di7Fjx6JTp05qCxjyfH19kZycjHHjxiEqKgo1atTAlClTMGXKFNHPc+jQIaP3nQIfYtecnDmHW3BO3xo0hBDrYa0ZH03h4eFCUbKqunXr4ocfftD72JYtW+Lo0aOynzsmJkbnfWfPnpXVJgU+xDYwpazVn3teStU7uyl99kgAT9ZmIUSsw5HdKPtDDFLOlChXyp/5Wm4Hm0ZLlZOTgy+++AKffPIJ/vjjDxrqIkSqgJkrLN0Fo6NMj3bGDnJjUw8Kw2e0dACRw1jr+DiCgwcPYvPmzdi1axfCwsLwwgsv4NNP5f0uU+BD7IOBn3zErPBLiCYKeAgxndu3b2Pr1q3YvHkzCgoKMGTIEJSWlmLnzp2IiIiQ3a71batNiJH9u3w6/l0+XdS5Y7hwYTFEQoCKYJgCYmIKfI2PITd79dxzzyEiIgIXL17EunXrcPfuXWHzU0NRxofYPX73djGfztfknjVtZ4jF8MOa/M8Bd/7xkGCL7joeoY4//9DwhQBodhcxXJkS4Azand2InbEySUlJmDhxIt5++200atTIqG1T4ENsl4SNTIEnb1xM5BudreEzVRS86Xc4shsA6YEL/3MTm2qfPz+EWJNffvkFmzdvRlRUFJo2bYrXX38dL7/8slHapsCH2A5ddTwi63vsNeAh4mQtmQqAMjXEepQrGZxsYDq7JURHRyM6Ohpr1qzBV199hc2bN2PKlClQKpVITk5GaGgovL29ZbVNNT7ErhTn56A4P0frfZN82mCSTxvzdohYBYV3dQTMXCF5Fp/Unxf6GSNSUI1P1Tw9PTFixAgcO3YMaWlpmDp1KhYvXoyAgAD0799fVpsU+BC7oqjmC0U1X4Pbodk69oW+n4TYviZNmmDp0qW4ffs2vvzyS9nt0FAXsSt8tscYwY+todoe45N6Tdcer8go2f/ncGIMNNQlj7OzMwYOHIiBAwfKejwFPsQuFRfkAQAUXk/GgNedXA0AEDsRouTETgCAW/QLxuwasVH8EJZqMKRZME91ZEQKWsDQMijwIXbt9txRwnR21/YD1e6raviD3sSIqo3sBgD1n5uJHSsKpvlgSHXKO19Ebe+zCQmxNRT4ELtWa8ZKg9ug/bwcm5RVvWnKO5GiXMkMWsfHUYe6DEWBD7Ed/Ho9Eren4Hcr73kpVdbT2uN+XsQ4xNQAUaaH6MIYAzMgeNG2YzqpGgU+xHZoLlgoMgDqcuqwafpDHIrUBSL5dYMocCa6KJXMoDodqvGRh6azE9vDlAZvSgpI35A0uVmkkD0ijmdN7lmdQY/mfl786tCEEOtDGR9in1QDo8cZIkOnunf5/WeDu0XsE18HxM/8GtTIjzI9pEqMMYOGq2ioSx4KfAgBAI6r+Jf+kBAD8Ov4/LP3Rwv3hNgCpjSwxoeGumShoS5i//QNjXHck6BH29fEbmQtmSrU3ZgKa9GdipkJsXKU8SFEojFcuLCmC7FO2uq3/tx5FgAQuqDiPn6BSgpUiKVQcbNlUOBDHJauzUyrQltDWB8xheqayxnwK3KbYh+vgJkraFYXqZKh8zSMMMfDIdFQFyHEpkmdnWfsxxNCbAtlfIjDkZvpUWtDwmq+xDSMee35toyZ/dHM9JjiOYhto1ldlkGBDyHE7mgGF1UFSVIXJ5SCn+LO14UpvKtT8EMAUI2PpVDgQwixG3IDCm0bkAKmCVIo6CHEsijwIYTYFM0Mij5Sh8M0zzdGkKIti8QXPvMb31Iw5JhoHR/LoMCHEGLzdAUOhtZiqT7O0OBE9fF8/c+Bxk8BAGJTDxrUNrFRBgY+oMBHFgp8CCE2TUxAUpz30CoL0SngcWxKxsAZUKCspOJmWSjwIUQbjqPtK6yU3MUj+QCJH2ZauPBnSe3xQ2ymXMeJCp8JMT0KfAgxAL/6L78YniWZ443ZksTU34iZMs4PM62ZKe35dRVAG4OUuiViPxgzsMaHPpzJQoEPIXbClG/M1oB/Xdz5ikwNLLTVxOHIbgBMM0ylOiRnr99H8gQVN1sGBT6ESMSd/1nY34n2eZJO841d6tYO2q655iwpTarfM01S6396XzktPM5Y+F3dGZ5kES0V2BFi7yjwIUQXfpd2jXSyW/QL9GlcJtUAg/+/OaZ0VxWg8s+tmc0xdUG0tqCP7yvt9WX/lEqAM2gBQyN2xoFQ4EOIRBT0iCclcMhaMlX2mzz/OEO/N1KGr4xRiKzv9VLAY/9oywrLoMCHEGJ01jh13BTk1uNQNocQy6Hd2QkhRiU36NFVn2MLFN7VRb9uhXd1hC741KZfLzEOpjT8Zkp79+5F+/bt4eHhAX9/fwwePFjt/lu3buH555+Hl5cX/P39MXHiRJSUlKidk5aWhpiYGHh4eKB27dqYP3++xTNVlPEhhBiFPUzJNnRGlb0vKUCMS6lkBtb4mC6A2LlzJ0aNGoXExER069YNjDGkpaUJ95eXl6Nv376oVasWjh07hvv372P48OFgjGHdunUAgNzcXPTs2RNdu3ZFSkoKrly5gri4OHh5eWHq1Kkm63tVKPAhxM6Ycrq1KaXPHmnxoZ/02SMBqNcMGXvYjmrECM9ap7OXlZVh0qRJWLZsGUaOHCkcb9KkifD/pKQkXLx4Eenp6QgJCQEArFixAnFxcVi4cCF8fHzwxRdfoKioCFu3boVCoUCLFi1w5coVrFy5ElOmTAHHTyAxMxrqIsQA3Pmfn6wr4+A2shsGZXssHfTwfQiYuUIYupL6veWvAR98EmIOubm5arfi4mKD2jt9+jTu3LkDJycntG3bFsHBwejTpw8uXLggnHPixAm0aNFCCHoAoHfv3iguLkZqaqpwTkxMDBQKhdo5d+/exY0bNwzqoyEo8CHEAKxFd6tby6f3ldPCWjO2xJqCheK8h0JmRvX/YvW+clpS3Q9xTHzGx5AbAISGhsLX11e4LVq0yKB+/f333wCAhIQEzJ49Gz/88ANq1KiBmJgYPHjwAACQmZmJwMBAtcfVqFEDbm5uyMzM1HkO/zV/jiVQ4EMIIYRYgJIxg28AkJ6ejpycHOE2a9Ysrc+XkJAAjuP03k6dOgXl4wWC3nvvPbzwwguIjIzEli1bwHEcvvnmG6E9bUNVjDG145rn8IXNlhrmAiwc+CxatAjt2rWDt7c3AgICMHDgQFy+fFntHMYYEhISEBISAg8PD8TGxqql2wgxOY57spihBvpUb7gDjZ/CgcZPWWWWamLHJwWYfOYnffZIoRaIEGvg4+OjdlMdWlI1fvx4XLp0Se+tRYsWCA4OBgBEREQIj1UoFKhfvz5u3boFAAgKCqqUtcnOzkZpaamQ1dF2TlZWFgBUygSZk0UDnyNHjmDcuHE4efIkkpOTUVZWhl69eqGgoEA4Z+nSpVi5ciU+/PBDpKSkICgoCD179kReXp4Fe06I9ctaMlVYL8aaxaYeRGzqQasMJowxO4vqwIguxhrqEsvf3x9NmzbVe3N3d0dkZCQUCoVaIqK0tBQ3btxAWFgYACA6Ohrnz59HRkaGcE5SUhIUCgUiIyOFc44ePao2xT0pKQkhISEIDw834MoZhmOWnlCv4t9//0VAQACOHDmCLl26gDGGkJAQxMfHY+bMiq2Ui4uLERgYiCVLlmD06NFVtpmbmwtfX19k3bkFHx8fU78EYs+0/KpIyfYkN4sU/t/zUqoxeiSKuWYRyc182dosJ0d5nY4qNzcXAXXCkJOTY7L3DP59qcF/voCzm6fsdspLHuGvTcNM0tf4+Hh8++232Lx5M8LCwrBs2TL873//w59//okaNWqgvLwcbdq0QWBgIJYtW4YHDx4gLi4OAwcOFKaz5+TkoEmTJujWrRveffddXL16FXFxcZgzZw5NZ+fl5OQAAPz8/AAA169fR2ZmJnr16iWco1AoEBMTg+PHj2sNfIqLi9Uq2nNzc03ca+IwtOzdxW8o6Rb9QpUPN2ewQ4xPyNqI+F6LRSs4E2u1bNkyuLi44PXXX0dhYSHat2+PgwcPokaNGgAAZ2dn7N27F2PHjkWnTp3g4eGBoUOHYvny5UIbvr6+SE5Oxrhx4xAVFYUaNWpgypQpmDJliqVeFgArCnwYY5gyZQo6d+6MFi1aAHhS9a2tKvzmzZta21m0aBHmzZtn2s4SQhyOmOBWG32ZHnNs0EqsF1MygxYhNNU6PgDg6uqK5cuXqwUymurWrYsffvhBbzstW7bE0aNHjd09g1jNrK7x48fj3Llz+PLLLyvdp60qXFdF+KxZs9Sq29PT003SX+JAGNM6zAVY53R2TVRfYhrps0cKBc/aprzrmwavWRSvWgfE38fXaFEBvf3iNyk15Eaks4qMz4QJE7Bnzx4cPXoUderUEY4HBQUBqMj88FXmQEVVuK6KcIVCobOinRBJ7OSPilv0C5RR0EPsLuua52gbnhLTDnf+50rDZdqySZp7eRljN3hCiIUzPowxjB8/Hrt27cLBgwdRr149tfvr1auHoKAgJCcnC8dKSkpw5MgRdOzY0dzdJY6Az+5oC3r0ZH4c2eHIbla1+KBU5g4mDMkQUvbHvph7VhepYNGMz7hx47Bjxw58//338Pb2Fmp6fH194eHhAY7jEB8fj8TERDRq1AiNGjVCYmIiPD09MXToUEt2ndgbAwMafio27bhNxOCDLQpiHJtSyQAr3aTUnlk08NmwYQMAIDY2Vu34li1bEBcXBwCYMWMGCgsLMXbsWGRnZ6N9+/ZISkqCt7e3mXtLiG40K8f6GbrzuikcaPwUAEhevJGvBzo0fCEA29uQllRgynIwZblBjyfSWTTwEVOYxXEcEhISkJCQYPoOEUIks8YVl7Xhgwx7CBL4mqDeVyr+taZgjhBrZxXFzYTYi+K8hzR8IdIYLhyAcVZHFsMeAh5NFPDYNsr4WAYFPoQYmTXWb5himMeaXp+tMVYxuOr3gIIg28OUSgMDH6URe+M4KPAhhEimbUo2EY8fHrTG/ckIsXcU+BC79++yaQCAWtN1r0BqCtac+ZFb7yK8Fgp6DGKK7Axf8GzogpqTfNoAMN8QpCNj5eVg5QZkfAx4rCOzqk1KTYE2KSWimPDXwJoCH01i34A1X4MxapkcZWiGDyQ2shsAKmqbNIMKzXMM4SjX1VTMuUlp8Evr4OTqIbsdZWkhMr6ZYNK+2iPK+BACaN2A1BFUtRqwruDGEYe65NZJaQYzql/zbRkj4NEla8lUnetLaXstlPEh9o4CH2LVjrarKALtklLFkAwfsGgGMBynP5jRseebMUnZwd0aiMnkGPJabCkjMcmnTaWghL8+6bNHGrx+k6FZM9WFM/mfM2gMdelbVHOSTxshwOH7QgGP+dCsLsugwIeYlmZAIlGVAY+u51OlLRjivzawf6K6ZOFNTMdw4aIyCtY8JGeNQhd8KiqI46fty83qqAY3uvYLO9D4KZ31WumzR+oNfvgZZvY43d/aUeBjGRT4ENtX1fCU5v269uEyMUut8VPVGy4FPLpVlf0QM/yl2YZqjY+2x2t+P/jgxhSZso3shlDoToijoMCHmJYZhpL0srKaHWuc6WVu1rh1hD7G+p7x7agGQsa6BnKzNcV5D0U9NmvJVAC0NYuxUcbHMijwIbbHyoIZOTTf8Bw5ELJVVQ0hGcpYW2wEzFxhcIDFBzwUABkXLWBoGRT4ENthBwEPsZ1MjybNzI8xAgp9zF1zIyYTRwGPcSmV5YABgY+SMj6yUOBDTEs1WLH0sBchRmDqwE3X7CxTs9WAlBCpKPAhxArQ5qaEZ+lZgMR8qMbHMijwIcRKqE5bJsQa2VphurWjwMcyKPAhxEpQ/QSxdsYquCbEkijwIeYnZ9FAByps5t9c+B28CbEWfMBjrA1RHV55OZiTAVkb2qRUFgp8iPk4UPBiCP7NxdAVfwkxFQp4jIMxw2Z1MUaBjxwU+NgDudsuGHu7BsZo5pYR0Z5JhBBifBT42ANde1FZgrGzOpQlohlfhNgpplQalvGhBQxlocDHHlBwYPdoqwtC7A8zcAFDmtUlDwU+ts4agh6xfdDMTInZSZ2ooQCIWBsqdCa2xsnSHSAiMaYeDGh+bQv4PnOc7uE4W3xdFlCc99Am11IRViXWQeFdnYI64jAq9uoy7Eako4yPrRETFIjdJsJYtUBiszRU+EyqYIvBnKNzi35B7Wv6HopHQ12WQRkfR2ENmRRr6IOdMTTzU5z3EGO4cGHqvKnpGw7hd/4mtoX/+bHVLKQl8Ss3G3Ij0lHGxxppzs4yJFjQ91hzzwKr6vVwHAVGMkmd+aX6BmWOafP8dhz6VqdWvU/b1gh8YEQrXFunw5HdANCqzsT6UeBjzaw9CJDaP2t/PTZOTOGzvk/kpiyc5vcfKzmxU8j6TPJpA0B/4KXaF809zCi7YB347x8f+BDxlMpycDTUZXYU+FgjbZkRR5jtZM+vzYwMDWD4AmTN2g1jUG2TVqS2L5TpkY6VKwHOgMCnnIqb5aDAxxzkDimpBjvGCgqMMbxFAYpN4AMg7vzPkqYa28q0ZNop3DrQdHZia6i42dL4oEZfMEGzoYgBLPmGlNwsEsnNIi32/MR8Dkd2o+EuiRgzsLjZhHt1XblyBQMGDIC/vz98fHzQqVMnHDp0SO2cW7du4fnnn4eXlxf8/f0xceJElJSUqJ2TlpaGmJgYeHh4oHbt2pg/fz6YhT88U+BjSrrW3hET7JgKv4YOZW1IFfiCZEP0vJSKnpdSjdAb3fhaIWIZrEV3yvbIZM2zuvr27YuysjIcPHgQqampaNOmDfr164fMzEwAQHl5Ofr27YuCggIcO3YMX331FXbu3ImpU5/MzszNzUXPnj0REhKClJQUrFu3DsuXL8fKlStN1m8xaKiLEEIIIYJ79+7h2rVr2Lx5M1q1agUAWLx4MdavX48LFy4gKCgISUlJuHjxItLT0xESEgIAWLFiBeLi4rBw4UL4+Pjgiy++QFFREbZu3QqFQoEWLVrgypUrWLlyJaZMmQLOQqMZlPGxJrqyQMbIztBqyUQiW5k2TkXSxFYZK+OTm5urdisuLjaoXzVr1kSzZs2wfft2FBQUoKysDB999BECAwMRGVkxdH3ixAm0aNFCCHoAoHfv3iguLkZqaqpwTkxMDBQKhdo5d+/exY0bNwzqoyEo42MttA0/GTMYkVrUrGstIao3sl/s8QwR7snnoQONnwIA9L5y2hI9IjaEZnVJx5Tlhs3qehz4hIaGqh2fO3cuEhISZLfLcRySk5MxYMAAeHt7w8nJCYGBgdi/fz+qV68OAMjMzERgYKDa42rUqAE3NzdhOCwzMxPh4eFq5/CPyczMRL169WT30RB2H/jwRVR5eXmWeHLx51pD3Y1qUKMv8LF0P4lxMY0psSqBT0H5k0+U1q2izzTDixiKf68wSwFueSkMepbyUgBAeno6fHx8hMOqGRZVCQkJmDdvnt4mU1JSEBkZibFjxyIgIAC//PILPDw88Mknn6Bfv35ISUlBcHAwAGgdqmKMqR3XPIe/rpYa5gIcIPDhf4gbNG1u4Z4QYsN8fS3dA0LMKi8vD74m+rl3c3NDUFAQMi/+1+C2goKC4O/vD3d39yrPHT9+PF555RW954SHh+PgwYP44YcfkJ2dLQRU69evR3JyMrZt24Z33nkHQUFB+O2339Qem52djdLSUiGrExQUJGR/eFlZWQBQKVtkTnYf+ISEhODixYuIiIioFBUT48jNzUVoaChdXxOh62s6dG1NyxavL2MMeXl5arUrxubu7o7r169Xmvoth5ubm6igBwD8/f3h7+9f5XmPHj0CADg5qZcBOzk5Qfl4R/jo6GgsXLgQGRkZQgYoKSkJCoVCqAOKjo7Gu+++i5KSEri5uQnnhISEVBoCMyeOWXpCvRnk5ubC19cXOTk5NvPLZ0vo+poWXV/ToWtrWnR9bdO9e/fQtGlTxMTEYM6cOfDw8MDHH3+MNWvWICUlBa1bt0Z5eTnatGmDwMBALFu2DA8ePEBcXBwGDhyIdevWAQBycnLQpEkTdOvWDe+++y6uXr2KuLg4zJkzR23au7nRrC5CCCGECPz9/bF//37k5+ejW7duiIqKwrFjx/D999+jdevWAABnZ2fs3bsX7u7u6NSpE4YMGYKBAwdi+fLlQju+vr5ITk7G7du3ERUVhbFjx2LKlCmYMmWKpV4aAMr4ECOg62tadH1Nh66tadH1JdbIITI+CoUCc+fO1VnpTgxD19e06PqaDl1b06LrS6yRQ2R8CCGEEEIAB8n4EEIIIYQAFPgQQgghxIFQ4EMIIYQQh0GBDyGEEEIchs0GPosWLUK7du3g7e2NgIAADBw4EJcvXxbuLy0txcyZM9GyZUt4eXkhJCQEb7zxBu7evavWTnFxMSZMmAB/f394eXmhf//+uH37trlfjtWp6vpqGj16NDiOw+rVq9WO0/WtTOy1vXTpEvr37w9fX194e3ujQ4cOuHXrlnA/XVvtxFzf/Px8jB8/HnXq1IGHhweaNWuGDRs2qJ1D11e7DRs2oFWrVvDx8YGPjw+io6Oxb98+4X7GGBISEhASEgIPDw/ExsbiwoULam3QtSUWxWxU79692ZYtW9j58+fZ2bNnWd++fVndunVZfn4+Y4yxhw8fsh49erCvv/6a/fnnn+zEiROsffv2LDIyUq2dMWPGsNq1a7Pk5GR2+vRp1rVrV9a6dWtWVlZmiZdlNaq6vqp2797NWrduzUJCQtiqVavU7qPrW5mYa3vt2jXm5+fHpk+fzk6fPs3++usv9sMPP7B//vlHOIeurXZiru9bb73FGjRowA4dOsSuX7/OPvroI+bs7My+++474Ry6vtrt2bOH7d27l12+fJldvnyZvfvuu8zV1ZWdP3+eMcbY4sWLmbe3N9u5cydLS0tjL7/8MgsODma5ublCG3RtiSXZbOCjKSsriwFgR44c0XnO77//zgCwmzdvMsYqgiNXV1f21VdfCefcuXOHOTk5sf3795u8z7ZE1/W9ffs2q127Njt//jwLCwtTC3zo+oqj7dq+/PLL7LXXXtP5GLq24mm7vs2bN2fz589XO++pp55is2fPZozR9ZWqRo0a7JNPPmFKpZIFBQWxxYsXC/cVFRUxX19ftnHjRsYYXVtieTY71KUpJycHAODn56f3HI7jUL16dQBAamoqSktL0atXL+GckJAQtGjRAsePHzdpf22NtuurVCrx+uuvY/r06WjevHmlx9D1FUfz2iqVSuzduxeNGzdG7969ERAQgPbt2+O7774THkPXVjxtP7udO3fGnj17cOfOHTDGcOjQIVy5cgW9e/cGQNdXrPLycnz11VcoKChAdHQ0rl+/jszMTLXrplAoEBMTI1w3urbE0uwi8GGMYcqUKejcuTNatGih9ZyioiK88847GDp0qLB0emZmJtzc3FCjRg21cwMDA5GZmWnyftsKXdd3yZIlcHFxwcSJE7U+jq5v1bRd26ysLOTn52Px4sV49tlnkZSUhEGDBmHw4ME4cuQIALq2Yun62V27di0iIiJQp04duLm54dlnn8X69evRuXNnAHR9q5KWloZq1apBoVBgzJgx2L17NyIiIoRrExgYqHa+6nWja0sszcXSHTCG8ePH49y5czh27JjW+0tLS/HKK69AqVRi/fr1VbbHGAPHccbups3Sdn1TU1OxZs0anD59WvK1ouv7hLZrq1QqAQADBgzA5MmTAQBt2rTB8ePHsXHjRsTExOhsj66tOl1/G9auXYuTJ09iz549CAsLw9GjRzF27FgEBwejR48eOtuj61uhSZMmOHv2LB4+fIidO3di+PDhQlAOoNI1EnPd6NoSc7H5jM+ECROwZ88eHDp0CHXq1Kl0f2lpKYYMGYLr168jOTlZbaO8oKAglJSUIDs7W+0xWVlZlT6xOCpd1/eXX35BVlYW6tatCxcXF7i4uODmzZuYOnUqwsPDAdD1rYqua+vv7w8XFxdERESond+sWTNhVhdd26rpur6FhYV49913sXLlSjz//PNo1aoVxo8fj5dfflnYWZqur35ubm5o2LAhoqKisGjRIrRu3Rpr1qxBUFAQAFTK3KheN7q2xNJsNvBhjGH8+PHYtWsXDh48iHr16lU6hw96rl69ip9++gk1a9ZUuz8yMhKurq5ITk4WjmVkZOD8+fPo2LGjyV+DNavq+r7++us4d+4czp49K9xCQkIwffp0HDhwAABdX12qurZubm5o165dpSnYV65cQVhYGAC6tvpUdX1LS0tRWloKJyf1P3/Ozs5Cto2urzSMMRQXF6NevXoICgpSu24lJSU4cuSIcN3o2hKLM389tXG8/fbbzNfXlx0+fJhlZGQIt0ePHjHGGCstLWX9+/dnderUYWfPnlU7p7i4WGhnzJgxrE6dOuynn35ip0+fZt26daNplazq66uN5qwuxuj6aiPm2u7atYu5urqyTZs2satXr7J169YxZ2dn9ssvvwjn0LXVTsz1jYmJYc2bN2eHDh1if//9N9uyZQtzd3dn69evF86h66vdrFmz2NGjR9n169fZuXPn2LvvvsucnJxYUlISY6xiOruvry/btWsXS0tLY6+++qrW6ex0bYml2GzgA0DrbcuWLYwxxq5fv67znEOHDgntFBYWsvHjxzM/Pz/m4eHB+vXrx27dumWZF2VFqrq+2mgLfOj6Vib22n766aesYcOGzN3dnbVu3VptjRnG6NrqIub6ZmRksLi4OBYSEsLc3d1ZkyZN2IoVK5hSqRTOoeur3YgRI1hYWBhzc3NjtWrVYt27dxeCHsYYUyqVbO7cuSwoKIgpFArWpUsXlpaWptYGXVtiSRxjjJkjs0QIIYQQYmk2W+NDCCGEECIVBT6EEEIIcRgU+BBCCCHEYVDgQwghhBCHQYEPIYQQQhwGBT6EEEIIcRgU+BBCCCHEYVDgQ2xCbGws4uPj7ep54+LiMHDgQIPaCA8PB8dx4DgODx8+1Hne1q1bUb16dYOei+gWFxcnfB++++47S3eHEKIHBT6E6LFr1y588MEHwtfh4eFYvXq15Tqkxfz585GRkQFfX19Ld8XuHT58WGuQuWbNGmRkZFimU4QQSVws3QFCrJmfn5+lu1Alb29vYVdsSystLYWrq6ulu2F2vr6+FHgSYiMo40NsUnZ2Nt544w3UqFEDnp6e6NOnD65evSrczw/tHDhwAM2aNUO1atXw7LPPqn0qLysrw8SJE1G9enXUrFkTM2fOxPDhw9WGn1SHumJjY3Hz5k1MnjxZGNYAgISEBLRp00atf6tXr0Z4eLjwdXl5OaZMmSI814wZM6C5WwxjDEuXLkX9+vXh4eGB1q1b49tvv5V1fbZu3Yq6devC09MTgwYNwv379yud87///Q+RkZFwd3dH/fr1MW/ePJSVlQn3//nnn+jcuTPc3d0RERGBn376SW0o58aNG+A4Dv/9738RGxsLd3d3fP755wCALVu2oFmzZnB3d0fTpk2xfv16tee+c+cOXn75ZdSoUQM1a9bEgAEDcOPGDeH+w4cP4+mnn4aXlxeqV6+OTp064ebNm6Jee1Wva+XKlWjZsiW8vLwQGhqKsWPHIj8/X7j/5s2beP7551GjRg14eXmhefPm+PHHH3Hjxg107doVAFCjRg1wHIe4uDhRfSKEWA8KfIhNiouLw6lTp7Bnzx6cOHECjDE899xzKC0tFc559OgRli9fjs8++wxHjx7FrVu3MG3aNOH+JUuW4IsvvsCWLVvw66+/Ijc3V299xq5du1CnTh1haEnK0MaKFSuwefNmfPrppzh27BgePHiA3bt3q50ze/ZsbNmyBRs2bMCFCxcwefJkvPbaazhy5Ij4CwPgt99+w4gRIzB27FicPXsWXbt2xYIFC9TOOXDgAF577TVMnDgRFy9exEcffYStW7di4cKFAAClUomBAwfC09MTv/32GzZt2oT33ntP6/PNnDkTEydOxKVLl9C7d298/PHHeO+997Bw4UJcunQJiYmJeP/997Ft2zYAFd+Xrl27olq1ajh69CiOHTsmBKYlJSUoKyvDwIEDERMTg3PnzuHEiRP4z3/+IwSa+lT1ugDAyckJa9euxfnz57Ft2zYcPHgQM2bMEO4fN24ciouLcfToUaSlpWHJkiWoVq0aQkNDsXPnTgDA5cuXkZGRgTVr1kj63hBCrIBFt0glRKSYmBg2adIkxhhjV65cYQDYr7/+Ktx/79495uHhwf773/8yxhjbsmULA8CuXbsmnPN///d/LDAwUPg6MDCQLVu2TPi6rKyM1a1blw0YMEDr8zKmfQf6uXPnstatW6sdW7VqFQsLCxO+Dg4OZosXLxa+Li0tZXXq1BGeKz8/n7m7u7Pjx4+rtTNy5Ej26quv6rwu2vrz6quvsmeffVbt2Msvv8x8fX2Fr5955hmWmJiods5nn33GgoODGWOM7du3j7m4uLCMjAzh/uTkZAaA7d69mzHG2PXr1xkAtnr1arV2QkND2Y4dO9SOffDBByw6OpoxVrHrfJMmTdR2Qi8uLmYeHh7swIED7P79+wwAO3z4sM7XrUtVr0ub//73v6xmzZrC1y1btmQJCQlazz106BADwLKzs7Xer3p9CCHWiWp8iM25dOkSXFxc0L59e+FYzZo10aRJE1y6dEk45unpiQYNGghfBwcHIysrCwCQk5ODf/75B08//bRwv7OzMyIjI6FUKo3a35ycHGRkZCA6Olo45uLigqioKGG46+LFiygqKkLPnj3VHltSUoK2bdtKer5Lly5h0KBBaseio6Oxf/9+4evU1FSkpKSoZULKy8tRVFSER48e4fLlywgNDVWrHVK9VqqioqKE///7779IT0/HyJEjMWrUKOF4WVmZUAOTmpqKa9euwdvbW62doqIi/PXXX+jVqxfi4uLQu3dv9OzZEz169MCQIUMQHBxc5Wuv6nV5enri0KFDSExMxMWLF5Gbm4uysjIUFRWhoKAAXl5emDhxIt5++20kJSWhR48eeOGFF9CqVasqn5sQYhso8CE2h2nUxqgeVx0O0Syy5Tiu0mM1h090ta2Pk5NTpcepDrmJwQdbe/fuRe3atdXuUygUktoS8xqUSiXmzZuHwYMHV7rP3d290rXUx8vLS61dAPj444/VAlOgIrDkz4mMjMQXX3xRqa1atWoBqKgRmjhxIvbv34+vv/4as2fPRnJyMjp06GDQ67p58yaee+45jBkzBh988AH8/Pxw7NgxjBw5UvievfXWW+jduzf27t2LpKQkLFq0CCtWrMCECRNEXQ9CiHWjwIfYnIiICJSVleG3335Dx44dAQD379/HlStX0KxZM1Ft+Pr6IjAwEL///jueeeYZABWZgTNnzlQqVFbl5uaG8vJytWO1atVCZmamWrBw9uxZtecKDg7GyZMn0aVLFwAVGZDU1FQ89dRTwmtSKBS4desWYmJiRL0GXSIiInDy5Em1Y5pfP/XUU7h8+TIaNmyotY2mTZvi1q1b+OeffxAYGAgASElJqfK5AwMDUbt2bfz9998YNmyY1nOeeuopfP311wgICICPj4/Ottq2bYu2bdti1qxZiI6Oxo4dO6oMfKp6XadOnUJZWRlWrFgBJ6eKEsf//ve/lc4LDQ3FmDFjMGbMGMyaNQsff/wxJkyYADc3NwCo9DNACLEdFPgQm9OoUSMMGDAAo0aNwkcffQRvb2+88847qF27NgYMGCC6nQkTJmDRokVo2LAhmjZtinXr1iE7O1tvpiM8PBxHjx7FK6+8AoVCAX9/f8TGxuLff//F0qVL8eKLL2L//v3Yt2+f2pv6pEmTsHjxYjRq1AjNmjXDypUr1daC8fb2xrRp0zB58mQolUp07twZubm5OH78OKpVq4bhw4eLfl0TJ05Ex44dsXTpUgwcOBBJSUlqw1wAMGfOHPTr1w+hoaF46aWX4OTkhHPnziEtLQ0LFixAz5490aBBAwwfPhxLly5FXl6eUNxcVSYoISEBEydOhI+PD/r06YPi4mKcOnUK2dnZmDJlCoYNG4Zly5ZhwIABmD9/PurUqYNbt25h165dmD59OkpLS7Fp0yb0798fISEhuHz5Mq5cuYI33nijytde1etq0KABysrKsG7dOjz//PP49ddfsXHjRrU24uPj0adPHzRu3BjZ2dk4ePCgEFCHhYWB4zj88MMPeO655+Dh4YFq1aqJ/t4QQqyAxaqLCJFAs8j4wYMH7PXXX2e+vr7Mw8OD9e7dm125ckW4f8uWLWrFvIwxtnv3bqb6I19aWsrGjx/PfHx8WI0aNdjMmTPZSy+9xF555RWdz3vixAnWqlUrplAo1NrasGEDCw0NZV5eXuyNN95gCxcuVCtuLi0tZZMmTWI+Pj6sevXqbMqUKeyNN95QK6RWKpVszZo1rEmTJszV1ZXVqlWL9e7dmx05ckTnddFW3MxYRQFxnTp1mIeHB3v++efZ8uXLK12P/fv3s44dOzIPDw/m4+PDnn76abZp0ybh/kuXLrFOnToxNzc31rRpU/a///2PAWD79+9njD0pbj5z5kyl5//iiy9YmzZtmJubG6tRowbr0qUL27Vrl3B/RkYGe+ONN5i/vz9TKBSsfv36bNSoUSwnJ4dlZmaygQMHsuDgYObm5sbCwsLYnDlzWHl5uc7rIOV1rVy5kgUHBws/N9u3b1crWB4/fjxr0KABUygUrFatWuz1119n9+7dEx4/f/58FhQUxDiOY8OHD1d7blBxMyFWj2NMRlEDIXZIqVSiWbNmGDJkiNpqzdYsPDwc8fHxZtnO49dff0Xnzp1x7do1taJx8gTHcdi9e7fBW5EQQkyH1vEhDuvmzZv4+OOPceXKFaSlpeHtt9/G9evXMXToUEt3TZKZM2eiWrVqyMnJMWq7u3fvRnJyMm7cuIGffvoJ//nPf9CpUycKerQYM2YMDXkRYiMo40McVnp6Ol555RWcP38ejDG0aNECixcvFgqQbcHNmzeF2Uj169cXCnaNYfv27fjggw+Qnp4Of39/9OjRAytWrEDNmjWN9hxSNW/eXOcKzh999JHOgmpTy8rKQm5uLoCKZRNUZ7oRQqwLBT6EEJuhGuhpCgwMrLQ2ECGEaKLAhxBCCCEOg2p8CCGEEOIwKPAhhBBCiMOgwIcQQgghDoMCH0IIIYQ4DAp8CCGEEOIwKPAhhBBCiMOgwIcQQgghDoMCH0IIIYQ4jP8HPAdIt0sWL58AAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Inspect variables and dimensions\n", + "print(ds.variables)\n", + "print(ds.dims)\n", + "\n", + "# Take a quick look at SST values\n", + "sst = ds['ANAL_TEMP']\n", + "print(sst)\n", + "\n", + "# Optional: visualize one time step\n", + "sst.isel(TIME1=0).plot()\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "d638946c-1f42-4c32-8344-e5c352f08f1b", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.DataArray 'ANAL_TEMP' (TIME1: 107, LAT: 361, AX005: 226)>\n",
+       "[8729702 values with dtype=float32]\n",
+       "Coordinates:\n",
+       "  * AX005    (AX005) float64 220.0 220.4 220.8 221.2 ... 308.8 309.2 309.6 310.0\n",
+       "  * LAT      (LAT) float32 15.0 15.12 15.25 15.38 ... 59.62 59.75 59.88 60.0\n",
+       "  * TIME1    (TIME1) datetime64[ns] 2005-06-09 2005-06-11 ... 2006-01-06\n",
+       "Attributes:\n",
+       "    long_name:      Analysis Temperature\n",
+       "    units:          Deg. C\n",
+       "    long_name_mod:  regrid: 0.4 deg on X\n",
+       "    history:        From CLASS_Descriptor
" + ], + "text/plain": [ + "\n", + "[8729702 values with dtype=float32]\n", + "Coordinates:\n", + " * AX005 (AX005) float64 220.0 220.4 220.8 221.2 ... 308.8 309.2 309.6 310.0\n", + " * LAT (LAT) float32 15.0 15.12 15.25 15.38 ... 59.62 59.75 59.88 60.0\n", + " * TIME1 (TIME1) datetime64[ns] 2005-06-09 2005-06-11 ... 2006-01-06\n", + "Attributes:\n", + " long_name: Analysis Temperature\n", + " units: Deg. C\n", + " long_name_mod: regrid: 0.4 deg on X\n", + " history: From CLASS_Descriptor" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# show one variable\n", + "ds[\"ANAL_TEMP\"] # e.g. the analysis temperature" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "2ef37d28-436f-4263-81da-0ec9e94d46c6", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:     (AX005: 226, LAT: 361, TIME1: 4, bnds: 2)\n",
+       "Coordinates:\n",
+       "  * AX005       (AX005) float64 220.0 220.4 220.8 221.2 ... 309.2 309.6 310.0\n",
+       "  * LAT         (LAT) float32 15.0 15.12 15.25 15.38 ... 59.62 59.75 59.88 60.0\n",
+       "  * TIME1       (TIME1) datetime64[ns] 2005-08-24 2005-08-26 ... 2005-08-30\n",
+       "Dimensions without coordinates: bnds\n",
+       "Data variables:\n",
+       "    TIME1_bnds  (TIME1, bnds) datetime64[ns] 2005-08-23 ... 2005-08-31\n",
+       "    ANAL_TEMP   (TIME1, LAT, AX005) float32 ...\n",
+       "Attributes:\n",
+       "    history:      FERRET V6.96   30-Nov-25\n",
+       "    Conventions:  CF-1.6
" + ], + "text/plain": [ + "\n", + "Dimensions: (AX005: 226, LAT: 361, TIME1: 4, bnds: 2)\n", + "Coordinates:\n", + " * AX005 (AX005) float64 220.0 220.4 220.8 221.2 ... 309.2 309.6 310.0\n", + " * LAT (LAT) float32 15.0 15.12 15.25 15.38 ... 59.62 59.75 59.88 60.0\n", + " * TIME1 (TIME1) datetime64[ns] 2005-08-24 2005-08-26 ... 2005-08-30\n", + "Dimensions without coordinates: bnds\n", + "Data variables:\n", + " TIME1_bnds (TIME1, bnds) datetime64[ns] 2005-08-23 ... 2005-08-31\n", + " ANAL_TEMP (TIME1, LAT, AX005) float32 ...\n", + "Attributes:\n", + " history: FERRET V6.96 30-Nov-25\n", + " Conventions: CF-1.6" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# select a time period e.g. for one specific hurricane\n", + "ds.sel(TIME1=slice(\"2005-08-23T00:00:00.000000000\", \"2005-08-31T00:00:00.000000000\")) # for hurricane katrina" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "35d89cba-af6f-4d5e-80fb-f687ac94bc85", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5a360fef-b0a7-40a2-900d-c9f812ddba23", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "py39-cartopy [conda env:py39-cartopy]", + "language": "python", + "name": "conda-env-py39-cartopy-py" + }, + "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.9.18" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}