diff --git a/notebooks/Gallery.ipynb b/notebooks/Gallery.ipynb index 19877578..fb826ed8 100644 --- a/notebooks/Gallery.ipynb +++ b/notebooks/Gallery.ipynb @@ -28,7 +28,8 @@ "- [Accessing ERA5 Data](./tutorial/Accessing_ERA5_Data.ipynb) (working as at 29/3/2025)\n", "- [Introduction to Pipelines](./tutorial/Data_Pipelines.ipynb) (working as at 29/3/2025)\n", "- [End-to-end CNN Training Example](./tutorial/CNN_model_training.ipynb) (working as at 29/3/2025)\n", - "- [Working with Multiple Data Sources](./tutorial/MultipleSources.ipynb) (working as at 29/3/2025)" + "- [Working with Multiple Data Sources](./tutorial/MultipleSources.ipynb) (working as at 29/3/2025)\n", + "- [Working with Climate Data](./tutorial/Working_with_Climate_Data.ipynb) (working as at 23/4/2025)" ] }, { @@ -99,7 +100,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.7" + "version": "3.12.9" } }, "nbformat": 4, diff --git a/notebooks/tutorial/Working_with_Climate_Data.ipynb b/notebooks/tutorial/Working_with_Climate_Data.ipynb new file mode 100644 index 00000000..b7ac0275 --- /dev/null +++ b/notebooks/tutorial/Working_with_Climate_Data.ipynb @@ -0,0 +1,5491 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "bd231875-890b-48f7-86fa-d89867799dfc", + "metadata": {}, + "source": [ + "# Working with Climate Data" + ] + }, + { + "cell_type": "markdown", + "id": "21fd252f-ab46-468b-9679-ffee05ce8572", + "metadata": {}, + "source": [ + "## Introduction\n", + "\n", + "This tutorial shows, in detail, how to make a PyEarthTools pipeline for loading climate data. It is a precursor step to creating an ML model for bias correction. It goes into a lot of detail in order to document how to process new data sources and what steps are involved. The end result is a re-usable data loading pipeline, which can then be re-used for many projects, based on the steps and assumptions shown here.\n", + "\n", + "In general terms, an \"analysis\" refers to estimates of current conditions, \"reanalysis\" refers to estimates of historical conditions, \"nowcasting\" refers to predictions of conditions which are expected in the next 90 minutes, \"short term\" refers to around 6-48 hours lead time and \"weather\" refers to 6 hours to 10 days lead time. A variety of terms like \"multi-week\", \"sub-seasonal\" and \"seasonal\" refer to the period from weeks to months. Lead times longer than a few months may be referred to as \"seasonal\" and merge into the climate time scales. Time scales for climate modelling are typically decades into the future or longer. These terms are not accepted universally, and some people may refer to any of these time scales as climate modelling. Within this tutorial, climate data is meant to describe predictions at least a year into the future.\n", + "\n", + "A \"reforecast\" is done to generate \"what the forecast would have been\". It is subtly different to a \"reanalysis\", because it includes the lead time component as well as the estimate at a point in time.\n", + "\n", + "Climate data uses non-Gregorian calendars, which also involve the use of date/time libraries which may be unfamiliar to many users. \n", + "\n", + "This tutorial will demonstrate the loading of climate data which includes both reforecasting outputs and predictions of the future, and merging it with reanalysis data for the purposes of validation of the reforecast component of the dataset.\n", + "\n", + "This will involve loading a climate prediction run from CMIP5 and ERA5 reanalysis data.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "c54f1450-09bb-4209-93b2-85f0b409cfb6", + "metadata": {}, + "outputs": [], + "source": [ + "import pyearthtools.data as petdata\n", + "import pyearthtools.pipeline as petpipe\n", + "import site_archive_nci" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "7b5b550e-3926-4c2e-8222-0a73d05b638c", + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "\n", + "# This builds an accessor that can be indexed by time, filtered according to the specified parameters\n", + "# Multiple institutions, scenarios and models are unsupported but the intention is to support that in future\n", + "# Models should generally be included in a pipeline of operations rather than used directly, but we will\n", + "# explore some of the functionality of this object regardless\n", + "cmip5_model1 = petdata.archive.CMIP5(institutions='BCC', scenarios=['rcp60'], models=['bcc-csm1-1'], interval='mon', variables='tas')" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "f50daa56-8456-44d1-bfd9-3266478dc634", + "metadata": {}, + "outputs": [], + "source": [ + "# With an exact time, you need to pick a time actually in the dataset, for fuzzy selection see the next cell\n", + "\n", + "# Under-specifying the datetime will request all source data which matches Jan 2010\n", + "# In this case, the data is monthy, with a pseudo-day-of-month of the 16th of the month\n", + "# Note for later - longitude is indexed from 0 to 360\n", + "ds_cmip_2010 = cmip5_model1['2010-01'] # Query data along primary dimension\n", + "ds_cmip_2010" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "1c0c36ce-cf18-4cad-8d89-9744003b0400", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 862kB\n",
+       "Dimensions:    (time: 24, bnds: 2, lat: 64, lon: 128)\n",
+       "Coordinates:\n",
+       "  * time       (time) object 192B 2010-01-16 12:00:00 ... 2011-12-16 12:00:00\n",
+       "  * lat        (lat) float64 512B -87.86 -85.1 -82.31 ... 82.31 85.1 87.86\n",
+       "  * lon        (lon) float64 1kB 0.0 2.812 5.625 8.438 ... 351.6 354.4 357.2\n",
+       "    height     float64 8B 2.0\n",
+       "Dimensions without coordinates: bnds\n",
+       "Data variables:\n",
+       "    time_bnds  (time, bnds) object 384B dask.array<chunksize=(24, 2), meta=np.ndarray>\n",
+       "    lat_bnds   (time, lat, bnds) float64 25kB dask.array<chunksize=(24, 64, 2), meta=np.ndarray>\n",
+       "    lon_bnds   (time, lon, bnds) float64 49kB dask.array<chunksize=(24, 128, 2), meta=np.ndarray>\n",
+       "    tas        (time, lat, lon) float32 786kB dask.array<chunksize=(24, 64, 128), meta=np.ndarray>\n",
+       "Attributes: (12/24)\n",
+       "    institution:            Beijing Climate Center(BCC),China Meteorological ...\n",
+       "    institute_id:           BCC\n",
+       "    experiment_id:          rcp60\n",
+       "    source:                 bcc-csm1-1:atmosphere:  BCC_AGCM2.1 (T42L26); lan...\n",
+       "    model_id:               bcc-csm1-1\n",
+       "    forcing:                Nat Ant GHG SD Oz Sl SS Ds BC OC\n",
+       "    ...                     ...\n",
+       "    table_id:               Table Amon (11 April 2011) 1cfdc7322cf2f4a3261482...\n",
+       "    title:                  bcc-csm1-1 model output prepared for CMIP5 RCP6\n",
+       "    parent_experiment:      historical\n",
+       "    modeling_realm:         atmos\n",
+       "    realization:            1\n",
+       "    cmor_version:           2.5.6
" + ], + "text/plain": [ + " Size: 862kB\n", + "Dimensions: (time: 24, bnds: 2, lat: 64, lon: 128)\n", + "Coordinates:\n", + " * time (time) object 192B 2010-01-16 12:00:00 ... 2011-12-16 12:00:00\n", + " * lat (lat) float64 512B -87.86 -85.1 -82.31 ... 82.31 85.1 87.86\n", + " * lon (lon) float64 1kB 0.0 2.812 5.625 8.438 ... 351.6 354.4 357.2\n", + " height float64 8B 2.0\n", + "Dimensions without coordinates: bnds\n", + "Data variables:\n", + " time_bnds (time, bnds) object 384B dask.array\n", + " lat_bnds (time, lat, bnds) float64 25kB dask.array\n", + " lon_bnds (time, lon, bnds) float64 49kB dask.array\n", + " tas (time, lat, lon) float32 786kB dask.array\n", + "Attributes: (12/24)\n", + " institution: Beijing Climate Center(BCC),China Meteorological ...\n", + " institute_id: BCC\n", + " experiment_id: rcp60\n", + " source: bcc-csm1-1:atmosphere: BCC_AGCM2.1 (T42L26); lan...\n", + " model_id: bcc-csm1-1\n", + " forcing: Nat Ant GHG SD Oz Sl SS Ds BC OC\n", + " ... ...\n", + " table_id: Table Amon (11 April 2011) 1cfdc7322cf2f4a3261482...\n", + " title: bcc-csm1-1 model output prepared for CMIP5 RCP6\n", + " parent_experiment: historical\n", + " modeling_realm: atmos\n", + " realization: 1\n", + " cmor_version: 2.5.6" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Here we show how to request a two-year data period. This is not the normal method\n", + "# for data access for the end-to-end example, but can be used to simply retrieve data\n", + "two_years_of_cmip = cmip5_model1.series(start='2010-01-01', end='2012-01-01', interval = (1, 'month'))\n", + "two_years_of_cmip" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "db420f81-f64e-40a8-89fe-9adc1b26652c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 6GB\n",
+       "Dimensions:    (longitude: 1440, latitude: 721, time: 744)\n",
+       "Coordinates:\n",
+       "  * longitude  (longitude) float32 6kB -180.0 -179.8 -179.5 ... 179.5 179.8\n",
+       "  * latitude   (latitude) float32 3kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0\n",
+       "  * time       (time) datetime64[ns] 6kB 2010-01-01 ... 2010-01-31T23:00:00\n",
+       "Data variables:\n",
+       "    2t         (time, latitude, longitude) float64 6GB dask.array<chunksize=(93, 91, 180), meta=np.ndarray>\n",
+       "Attributes:\n",
+       "    Conventions:  CF-1.6\n",
+       "    history:      2020-09-28 12:42:00 UTC+1000 by era5_replication_tools-1.2....\n",
+       "    license:      Licence to use Copernicus Products: https://apps.ecmwf.int/...\n",
+       "    summary:      ERA5 is the fifth generation ECMWF atmospheric reanalysis o...\n",
+       "    title:        ERA5 single-levels reanalysis 2m_temperature 20100101-20100131
" + ], + "text/plain": [ + " Size: 6GB\n", + "Dimensions: (longitude: 1440, latitude: 721, time: 744)\n", + "Coordinates:\n", + " * longitude (longitude) float32 6kB -180.0 -179.8 -179.5 ... 179.5 179.8\n", + " * latitude (latitude) float32 3kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0\n", + " * time (time) datetime64[ns] 6kB 2010-01-01 ... 2010-01-31T23:00:00\n", + "Data variables:\n", + " 2t (time, latitude, longitude) float64 6GB dask.array\n", + "Attributes:\n", + " Conventions: CF-1.6\n", + " history: 2020-09-28 12:42:00 UTC+1000 by era5_replication_tools-1.2....\n", + " license: Licence to use Copernicus Products: https://apps.ecmwf.int/...\n", + " summary: ERA5 is the fifth generation ECMWF atmospheric reanalysis o...\n", + " title: ERA5 single-levels reanalysis 2m_temperature 20100101-20100131" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Here we create an ERA5 accessor. When we pass the same year and date in, this time around\n", + "# we get 744 time-steps, due to the increased time resolution of ERA5. In order to relate\n", + "# ERA5 to CMIP, we will need to deal with this temporal mismatch. They also have different\n", + "# spatial grids. We'll look at how to do that next.\n", + "era5 = petdata.archive.ERA5('2t',)\n", + "era5['2010-01'] # Uncomment this line to see the returned dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "747b8726-e767-4f30-8ff6-f92463105b17", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 8MB\n",
+       "Dimensions:    (longitude: 1440, latitude: 721, time: 1)\n",
+       "Coordinates:\n",
+       "  * longitude  (longitude) float32 6kB -180.0 -179.8 -179.5 ... 179.5 179.8\n",
+       "  * latitude   (latitude) float32 3kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0\n",
+       "  * time       (time) datetime64[ns] 8B 2010-01-01\n",
+       "Data variables:\n",
+       "    2t         (time, latitude, longitude) float64 8MB dask.array<chunksize=(1, 721, 1440), meta=np.ndarray>\n",
+       "Attributes:\n",
+       "    Conventions:  CF-1.6\n",
+       "    history:      2020-09-05 13:01:05 UTC+1000 by era5_replication_tools-1.0....\n",
+       "    license:      Licence to use Copernicus Products: https://apps.ecmwf.int/...\n",
+       "    summary:      ERA5 is the fifth generation ECMWF atmospheric reanalysis o...\n",
+       "    title:        ERA5 single-levels monthly-averaged 2m_temperature 20100101...
" + ], + "text/plain": [ + " Size: 8MB\n", + "Dimensions: (longitude: 1440, latitude: 721, time: 1)\n", + "Coordinates:\n", + " * longitude (longitude) float32 6kB -180.0 -179.8 -179.5 ... 179.5 179.8\n", + " * latitude (latitude) float32 3kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0\n", + " * time (time) datetime64[ns] 8B 2010-01-01\n", + "Data variables:\n", + " 2t (time, latitude, longitude) float64 8MB dask.array\n", + "Attributes:\n", + " Conventions: CF-1.6\n", + " history: 2020-09-05 13:01:05 UTC+1000 by era5_replication_tools-1.0....\n", + " license: Licence to use Copernicus Products: https://apps.ecmwf.int/...\n", + " summary: ERA5 is the fifth generation ECMWF atmospheric reanalysis o...\n", + " title: ERA5 single-levels monthly-averaged 2m_temperature 20100101..." + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# In this case, we are lucky, and ERA5 actually has a monthly-average product on disk already\n", + "# With this small change, we can ensure we're working with temporally matched data\n", + "# Note for later - longitude is indexed from -180 to + 180\n", + "era5 = petdata.archive.ERA5('2t', product='monthly-averaged')\n", + "era5['2010-01'] # Uncomment this line to see the returned dataset" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "47d54afd-8719-4c29-be0a-9a232a0970e2", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Pipeline\n",
+       "\tDescription                    `pyearthtools.pipeline` Data Pipeline\n",
+       "\n",
+       "\n",
+       "\tInitialisation                 \n",
+       "\t\t exceptions_to_ignore           None\n",
+       "\t\t iterator                       None\n",
+       "\t\t sampler                        None\n",
+       "\tSteps                          \n",
+       "\t\t _.CMIP5                        {'CMIP5': {'institutions': "'BCC'", 'interval': "'mon'", 'models': "['bcc-csm1-1']", 'scenarios': "['rcp60']", 'variables': "'tas'"}}\n",
+       "\t\t ERA5                           {'ERA5': {'level_value': 'None', 'product': "'monthly-averaged'", 'variables': "['2t']"}}
" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "

Graph

" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "%3\n", + "\n", + "\n", + "\n", + "CMIP5_11061510-b6e2-43ce-a60d-7a8bab926a70\n", + "\n", + "_.CMIP5\n", + "\n", + "\n", + "\n", + "ERA5_311760f1-0243-4221-8fb7-23ec1eb23994\n", + "\n", + "ERA5\n", + "\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "pipe = petpipe.Pipeline(\n", + " (cmip5_model1, era5),\n", + " )\n", + "pipe" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "45533105-1d57-453c-9545-3ea0c6925d04", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "( Size: 37kB\n", + " Dimensions: (time: 1, bnds: 2, latitude: 64, longitude: 128)\n", + " Coordinates:\n", + " * time (time) object 8B 2010-01-16 12:00:00\n", + " * latitude (latitude) float64 512B -87.86 -85.1 -82.31 ... 82.31 85.1 87.86\n", + " * longitude (longitude) float64 1kB 0.0 2.812 5.625 ... 351.6 354.4 357.2\n", + " height float64 8B 2.0\n", + " Dimensions without coordinates: bnds\n", + " Data variables:\n", + " time_bnds (time, bnds) object 16B dask.array\n", + " lat_bnds (time, latitude, bnds) float64 1kB dask.array\n", + " lon_bnds (time, longitude, bnds) float64 2kB dask.array\n", + " tas (time, latitude, longitude) float32 33kB dask.array\n", + " Attributes: (12/24)\n", + " institution: Beijing Climate Center(BCC),China Meteorological ...\n", + " institute_id: BCC\n", + " experiment_id: rcp60\n", + " source: bcc-csm1-1:atmosphere: BCC_AGCM2.1 (T42L26); lan...\n", + " model_id: bcc-csm1-1\n", + " forcing: Nat Ant GHG SD Oz Sl SS Ds BC OC\n", + " ... ...\n", + " table_id: Table Amon (11 April 2011) 1cfdc7322cf2f4a3261482...\n", + " title: bcc-csm1-1 model output prepared for CMIP5 RCP6\n", + " parent_experiment: historical\n", + " modeling_realm: atmos\n", + " realization: 1\n", + " cmor_version: 2.5.6,\n", + " Size: 8MB\n", + " Dimensions: (longitude: 1440, latitude: 721, time: 1)\n", + " Coordinates:\n", + " * longitude (longitude) float32 6kB -180.0 -179.8 -179.5 ... 179.5 179.8\n", + " * latitude (latitude) float32 3kB 90.0 89.75 89.5 ... -89.5 -89.75 -90.0\n", + " * time (time) datetime64[ns] 8B 2010-01-01\n", + " Data variables:\n", + " 2t (time, latitude, longitude) float64 8MB dask.array\n", + " Attributes:\n", + " Conventions: CF-1.6\n", + " history: 2020-09-05 13:01:05 UTC+1000 by era5_replication_tools-1.0....\n", + " license: Licence to use Copernicus Products: https://apps.ecmwf.int/...\n", + " summary: ERA5 is the fifth generation ECMWF atmospheric reanalysis o...\n", + " title: ERA5 single-levels monthly-averaged 2m_temperature 20100101...)" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Here we can draw matched time steps from CMIP5 and ERA5 \n", + "# There are still some problems to iron out, but things are getting closer\n", + "pipe['2010-01'] # Uncomment this line to see a simple use of the pipeline so far\n", + "\n", + "# The problems which remain are:\n", + "# - The climate data is in a 'noleap' calendar whereas ERA5 uses familiar date-time indexing\n", + "# - The variables are in a different order in the two objects (time,lat,lon,height) vs (lon,lat,time) assumed to be at the surface\n", + "# - ERA5 is still much higher spatial resolution than the CMIP data, so the lat/lon grids still need to be matched\n", + "# - The actual latitudes and longitudes may not match nicely and a bit of interpolation might be needed" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "a4273339-fb58-4989-8f35-0200b485abda", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Pipeline\n",
+       "\tDescription                    `pyearthtools.pipeline` Data Pipeline\n",
+       "\n",
+       "\n",
+       "\tInitialisation                 \n",
+       "\t\t exceptions_to_ignore           None\n",
+       "\t\t iterator                       None\n",
+       "\t\t sampler                        None\n",
+       "\tSteps                          \n",
+       "\t\t _.CMIP5                        {'CMIP5': {'institutions': "'BCC'", 'interval': "'mon'", 'models': "['bcc-csm1-1']", 'scenarios': "['rcp60']", 'variables': "'tas'"}}\n",
+       "\t\t _recode_calendar.RecodeCalendar {'RecodeCalendar': {}}\n",
+       "\t\t ERA5                           {'ERA5': {'level_value': 'None', 'product': "'monthly-averaged'", 'variables': "['2t']"}}\n",
+       "\t\t join.Merge                     {'Merge': {'merge_kwargs': 'None'}}
" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "

Graph

" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "%3\n", + "\n", + "\n", + "\n", + "CMIP5_7cf31d7e-6159-4030-b0aa-d7d4c117859a\n", + "\n", + "_.CMIP5\n", + "\n", + "\n", + "\n", + "RecodeCalendar_b17ebf99-38c0-4aad-92d0-e62ca47a9a7e\n", + "\n", + "_recode_calendar.RecodeCalendar\n", + "\n", + "\n", + "\n", + "CMIP5_7cf31d7e-6159-4030-b0aa-d7d4c117859a->RecodeCalendar_b17ebf99-38c0-4aad-92d0-e62ca47a9a7e\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "Merge_2ccca2bd-479d-43be-934f-e28af3c8bcee\n", + "\n", + "join.Merge\n", + "\n", + "\n", + "\n", + "RecodeCalendar_b17ebf99-38c0-4aad-92d0-e62ca47a9a7e->Merge_2ccca2bd-479d-43be-934f-e28af3c8bcee\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "ERA5_d19ba70f-dac5-4023-913d-7b98dcb996fc\n", + "\n", + "ERA5\n", + "\n", + "\n", + "\n", + "ERA5_d19ba70f-dac5-4023-913d-7b98dcb996fc->Merge_2ccca2bd-479d-43be-934f-e28af3c8bcee\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "# Before all tackling those issues, we're going to look a bit more at how PyEarthTools can be used to relate\n", + "# the data sources to produce input data for a machine learning model\n", + "# The basic goal is to define all the operations needed so that at the end of the day, it's easy\n", + "# to pull samples off an end-to-end pipeline, either by date/time (most common) or by another simple index (less common)\n", + "\n", + "# Here we start by producing a pipeline which partly tackles the calendar issue\n", + "# and merges the two data sources into a single in-memory object, in this case \n", + "# an xarray Dataset\n", + "\n", + "pipe = petpipe.Pipeline(\n", + " cmip5_model1,\n", + " (petpipe.operations.xarray.RecodeCalendar(), era5),\n", + " petpipe.operations.xarray.Merge()\n", + ")\n", + "pipe" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "e949afde-1135-496a-9a99-f1f66703d6c7", + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "# # Turn off debugging output, which reports the risks involved with recoding the calendar\n", + "# # and may produce some dask performance warnings\n", + "\n", + "sample = pipe['2010-01'] " + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "c9675055-ea17-483d-a416-bca4705ba42f", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 29MB\n",
+       "Dimensions:    (time: 2, latitude: 785, longitude: 1552, bnds: 2)\n",
+       "Coordinates:\n",
+       "  * time       (time) datetime64[ns] 16B 2010-01-01 2010-01-16T12:00:00\n",
+       "  * latitude   (latitude) float64 6kB -90.0 -89.75 -89.5 ... 89.5 89.75 90.0\n",
+       "  * longitude  (longitude) float64 12kB -180.0 -179.8 -179.5 ... 354.4 357.2\n",
+       "    height     float64 8B 2.0\n",
+       "Dimensions without coordinates: bnds\n",
+       "Data variables:\n",
+       "    time_bnds  (time, bnds) object 32B dask.array<chunksize=(1, 2), meta=np.ndarray>\n",
+       "    lat_bnds   (time, latitude, bnds) float64 25kB dask.array<chunksize=(1, 64, 2), meta=np.ndarray>\n",
+       "    lon_bnds   (time, longitude, bnds) float64 50kB dask.array<chunksize=(1, 128, 2), meta=np.ndarray>\n",
+       "    tas        (time, latitude, longitude) float32 10MB dask.array<chunksize=(1, 64, 128), meta=np.ndarray>\n",
+       "    2t         (time, latitude, longitude) float64 19MB dask.array<chunksize=(1, 721, 1552), meta=np.ndarray>\n",
+       "Attributes: (12/24)\n",
+       "    institution:            Beijing Climate Center(BCC),China Meteorological ...\n",
+       "    institute_id:           BCC\n",
+       "    experiment_id:          rcp60\n",
+       "    source:                 bcc-csm1-1:atmosphere:  BCC_AGCM2.1 (T42L26); lan...\n",
+       "    model_id:               bcc-csm1-1\n",
+       "    forcing:                Nat Ant GHG SD Oz Sl SS Ds BC OC\n",
+       "    ...                     ...\n",
+       "    table_id:               Table Amon (11 April 2011) 1cfdc7322cf2f4a3261482...\n",
+       "    title:                  bcc-csm1-1 model output prepared for CMIP5 RCP6\n",
+       "    parent_experiment:      historical\n",
+       "    modeling_realm:         atmos\n",
+       "    realization:            1\n",
+       "    cmor_version:           2.5.6
" + ], + "text/plain": [ + " Size: 29MB\n", + "Dimensions: (time: 2, latitude: 785, longitude: 1552, bnds: 2)\n", + "Coordinates:\n", + " * time (time) datetime64[ns] 16B 2010-01-01 2010-01-16T12:00:00\n", + " * latitude (latitude) float64 6kB -90.0 -89.75 -89.5 ... 89.5 89.75 90.0\n", + " * longitude (longitude) float64 12kB -180.0 -179.8 -179.5 ... 354.4 357.2\n", + " height float64 8B 2.0\n", + "Dimensions without coordinates: bnds\n", + "Data variables:\n", + " time_bnds (time, bnds) object 32B dask.array\n", + " lat_bnds (time, latitude, bnds) float64 25kB dask.array\n", + " lon_bnds (time, longitude, bnds) float64 50kB dask.array\n", + " tas (time, latitude, longitude) float32 10MB dask.array\n", + " 2t (time, latitude, longitude) float64 19MB dask.array\n", + "Attributes: (12/24)\n", + " institution: Beijing Climate Center(BCC),China Meteorological ...\n", + " institute_id: BCC\n", + " experiment_id: rcp60\n", + " source: bcc-csm1-1:atmosphere: BCC_AGCM2.1 (T42L26); lan...\n", + " model_id: bcc-csm1-1\n", + " forcing: Nat Ant GHG SD Oz Sl SS Ds BC OC\n", + " ... ...\n", + " table_id: Table Amon (11 April 2011) 1cfdc7322cf2f4a3261482...\n", + " title: bcc-csm1-1 model output prepared for CMIP5 RCP6\n", + " parent_experiment: historical\n", + " modeling_realm: atmos\n", + " realization: 1\n", + " cmor_version: 2.5.6" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sample" + ] + }, + { + "cell_type": "markdown", + "id": "b16a3736-3836-41df-b343-e06013c490cc", + "metadata": {}, + "source": [ + "This shows how a pipeline can be used to request all matching elements of each source according to month. The CMIP calendar has been recoded, but there is still a mismatch between the two different monthly values. The CMIP data in this example nominally has a day-of-month of the 16th of each month (e.g. 2010-01-16) whereas the ERA5 calendar chooses the start of the month (e.g. 2010-01-01). This introduces a lot of uncertainty to users. Does this mean the ERA5 averaging window is actually from the middle of one month to the middle of the next, centered on the first of the month? Or is it the rest of the month specified? A misinterpretation here could result in a temporal error of up to a fornight between the data sources. It is necessary to consult the data dictionary for the data source to determine what applies. Nuances like this are why it's a good idea for the community to encode standard loading pipelines into open source tools, so that newcomers (or people who might simply overlook the issue) can rely on standardised routines.\n", + "\n", + "In this case, both of these averages are actually for the same time period, being all of January 2010. We could just bake a translation into the data loading, but perhaps someone will come up with some other dataset or use case where the assumption isn't valid. As a result, we will make the translation explicit by adding another operation, just like the calendar recode step." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "6fca6d0a-99a0-4ebc-a643-d819c79c472f", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
Pipeline\n",
+       "\tDescription                    `pyearthtools.pipeline` Data Pipeline\n",
+       "\n",
+       "\n",
+       "\tInitialisation                 \n",
+       "\t\t exceptions_to_ignore           None\n",
+       "\t\t iterator                       None\n",
+       "\t\t sampler                        None\n",
+       "\tSteps                          \n",
+       "\t\t _.CMIP5                        {'CMIP5': {'institutions': "'BCC'", 'interval': "'mon'", 'models': "['bcc-csm1-1']", 'scenarios': "['rcp60']", 'variables': "'tas'"}}\n",
+       "\t\t _recode_calendar.RecodeCalendar {'RecodeCalendar': {}}\n",
+       "\t\t _align_dates.AlignDates        {'AlignDates': {'to': "'01'"}}\n",
+       "\t\t ERA5                           {'ERA5': {'level_value': 'None', 'product': "'monthly-averaged'", 'variables': "['2t']"}}\n",
+       "\t\t coordinates.StandardLongitude  {'StandardLongitude': {'longitude_name': "'longitude'", 'type': "'0-360'"}}\n",
+       "\t\t join.GeospatialTimeSeriesMerge {'GeospatialTimeSeriesMerge': {'interpolation_method': "'nearest'", 'merge_kwargs': 'None', 'reference_dataset': 'None', 'reference_index': '0', 'time_dimension': "'time'"}}
" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/html": [ + "

Graph

" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/svg+xml": [ + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "%3\n", + "\n", + "\n", + "\n", + "CMIP5_5c74717e-b1f5-4646-b3b8-061de949797c\n", + "\n", + "_.CMIP5\n", + "\n", + "\n", + "\n", + "RecodeCalendar_df01d847-5a30-40ea-90ff-1cd041fb87b4\n", + "\n", + "_recode_calendar.RecodeCalendar\n", + "\n", + "\n", + "\n", + "CMIP5_5c74717e-b1f5-4646-b3b8-061de949797c->RecodeCalendar_df01d847-5a30-40ea-90ff-1cd041fb87b4\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "AlignDates_3df195ee-0992-4845-bd1a-2acc14560525\n", + "\n", + "_align_dates.AlignDates\n", + "\n", + "\n", + "\n", + "RecodeCalendar_df01d847-5a30-40ea-90ff-1cd041fb87b4->AlignDates_3df195ee-0992-4845-bd1a-2acc14560525\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "GeospatialTimeSeriesMerge_e29b58eb-7bf7-4d8e-a234-f3d947713cd0\n", + "\n", + "join.GeospatialTimeSeriesMerge\n", + "\n", + "\n", + "\n", + "AlignDates_3df195ee-0992-4845-bd1a-2acc14560525->GeospatialTimeSeriesMerge_e29b58eb-7bf7-4d8e-a234-f3d947713cd0\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "ERA5_99763955-f39d-4472-abb9-e20e2031a3f3\n", + "\n", + "ERA5\n", + "\n", + "\n", + "\n", + "StandardLongitude_59909a17-709f-4f22-ad56-4746701fba0b\n", + "\n", + "coordinates.StandardLongitude\n", + "\n", + "\n", + "\n", + "ERA5_99763955-f39d-4472-abb9-e20e2031a3f3->StandardLongitude_59909a17-709f-4f22-ad56-4746701fba0b\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "StandardLongitude_59909a17-709f-4f22-ad56-4746701fba0b->GeospatialTimeSeriesMerge_e29b58eb-7bf7-4d8e-a234-f3d947713cd0\n", + "\n", + "\n", + "\n", + "\n", + "\n" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "pipe1 = petpipe.Pipeline(\n", + " cmip5_model1,\n", + " petpipe.operations.xarray.RecodeCalendar(),\n", + " petpipe.operations.xarray.AlignDates(to=\"start_of_month\"), \n", + ")\n", + "\n", + "pipe2 = petpipe.Pipeline(\n", + " era5,\n", + " petdata.transforms.coordinates.StandardLongitude(type=\"0-360\"), \n", + " \n", + ")\n", + "\n", + "pipe3 = petpipe.Pipeline(\n", + " (pipe1, pipe2),\n", + " petpipe.operations.xarray.join.GeospatialTimeSeriesMerge(reference_index=0)\n", + ")\n", + "\n", + "pipe3" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "fddd9dcb-6607-488b-811e-4505c0848a4e", + "metadata": {}, + "outputs": [], + "source": [ + "%%capture\n", + "# Turn off debugging output, which reports the risks involved with recoding the calendar\n", + "# and may produce some dask performance warnings\n", + "\n", + "sample = pipe3['2010-01']" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "4f288208-9a43-4c97-ac15-1f215d9c492c", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset> Size: 103kB\n",
+       "Dimensions:    (time: 1, latitude: 64, bnds: 2, longitude: 128)\n",
+       "Coordinates:\n",
+       "    height     float64 8B 2.0\n",
+       "  * latitude   (latitude) float64 512B -87.86 -85.1 -82.31 ... 82.31 85.1 87.86\n",
+       "  * longitude  (longitude) float64 1kB 0.0 2.812 5.625 ... 351.6 354.4 357.2\n",
+       "  * time       (time) datetime64[ns] 8B 2010-01-01\n",
+       "Dimensions without coordinates: bnds\n",
+       "Data variables:\n",
+       "    lat_bnds   (time, latitude, bnds) float64 1kB dask.array<chunksize=(1, 64, 2), meta=np.ndarray>\n",
+       "    lon_bnds   (time, longitude, bnds) float64 2kB dask.array<chunksize=(1, 128, 2), meta=np.ndarray>\n",
+       "    tas        (time, latitude, longitude) float32 33kB dask.array<chunksize=(1, 64, 128), meta=np.ndarray>\n",
+       "    time_bnds  (time, bnds) object 16B dask.array<chunksize=(1, 2), meta=np.ndarray>\n",
+       "    2t         (time, latitude, longitude) float64 66kB dask.array<chunksize=(1, 64, 128), meta=np.ndarray>\n",
+       "Attributes: (12/24)\n",
+       "    institution:            Beijing Climate Center(BCC),China Meteorological ...\n",
+       "    institute_id:           BCC\n",
+       "    experiment_id:          rcp60\n",
+       "    source:                 bcc-csm1-1:atmosphere:  BCC_AGCM2.1 (T42L26); lan...\n",
+       "    model_id:               bcc-csm1-1\n",
+       "    forcing:                Nat Ant GHG SD Oz Sl SS Ds BC OC\n",
+       "    ...                     ...\n",
+       "    table_id:               Table Amon (11 April 2011) 1cfdc7322cf2f4a3261482...\n",
+       "    title:                  bcc-csm1-1 model output prepared for CMIP5 RCP6\n",
+       "    parent_experiment:      historical\n",
+       "    modeling_realm:         atmos\n",
+       "    realization:            1\n",
+       "    cmor_version:           2.5.6
" + ], + "text/plain": [ + " Size: 103kB\n", + "Dimensions: (time: 1, latitude: 64, bnds: 2, longitude: 128)\n", + "Coordinates:\n", + " height float64 8B 2.0\n", + " * latitude (latitude) float64 512B -87.86 -85.1 -82.31 ... 82.31 85.1 87.86\n", + " * longitude (longitude) float64 1kB 0.0 2.812 5.625 ... 351.6 354.4 357.2\n", + " * time (time) datetime64[ns] 8B 2010-01-01\n", + "Dimensions without coordinates: bnds\n", + "Data variables:\n", + " lat_bnds (time, latitude, bnds) float64 1kB dask.array\n", + " lon_bnds (time, longitude, bnds) float64 2kB dask.array\n", + " tas (time, latitude, longitude) float32 33kB dask.array\n", + " time_bnds (time, bnds) object 16B dask.array\n", + " 2t (time, latitude, longitude) float64 66kB dask.array\n", + "Attributes: (12/24)\n", + " institution: Beijing Climate Center(BCC),China Meteorological ...\n", + " institute_id: BCC\n", + " experiment_id: rcp60\n", + " source: bcc-csm1-1:atmosphere: BCC_AGCM2.1 (T42L26); lan...\n", + " model_id: bcc-csm1-1\n", + " forcing: Nat Ant GHG SD Oz Sl SS Ds BC OC\n", + " ... ...\n", + " table_id: Table Amon (11 April 2011) 1cfdc7322cf2f4a3261482...\n", + " title: bcc-csm1-1 model output prepared for CMIP5 RCP6\n", + " parent_experiment: historical\n", + " modeling_realm: atmos\n", + " realization: 1\n", + " cmor_version: 2.5.6" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "sample" + ] + }, + { + "cell_type": "markdown", + "id": "2fefe5f7-c413-4b23-a5aa-5cf4bfca0ead", + "metadata": {}, + "source": [ + "Aligning the dates mean that the merged pipeline samples can use a common index. When the dates were aligned, the variable sorting issue also disappeared, because xarray could now see how to match the variables across the sources more effectively. It is possible to control the sort order of variables in an array, but it is not necessary here any longer.\n", + "\n", + "The final issue pertains to the grid resolution of the two datasets. Let's do some plotting to make it clearer.\n", + "\n", + "The problems which remain are:\n", + "- ERA5 is still much higher spatial resolution than the CMIP data, so the lat/lon grids still need to be matched\n", + "- The actual latitudes and longitudes may not match nicely and a bit of interpolation might be needed" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "887ffe1a-a58a-45e7-aef6-9737778118ae", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 15, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlEAAAHFCAYAAADSY6wWAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAwN5JREFUeJzsnXl4FEX+/989M5nJ5CQhkIQznMqpCIqAK6ByuIjX/tYDl1tXRUVERFldCSoghyweq4gH4NdbwQtXFlaFFUHlEJVDXM5EIAaSkDtz9u+PmK5PTaYmM8kkmcDn9Tx5npqe6qrq6upOTfXn3W9N13UdDMMwDMMwTEiYGrsBDMMwDMMwTRGeRDEMwzAMw9QCnkQxDMMwDMPUAp5EMQzDMAzD1AKeRDEMwzAMw9QCnkQxDMMwDMPUAp5EMQzDMAzD1AKeRDEMwzAMw9QCnkQxDMMwDMPUAp5EnSFkZmZC0zScOnUqbGVOmDABGRkZ9d6eefPm4cMPP6xVPXWhqKgIc+fOxZAhQ5CWloa4uDj06tULCxYsQEVFRdDlvP322zj//PMRHR2NVq1aYdq0aSgpKalxv5UrV0LTNOMvnOfuww8/lMrevn17yGUcP34cmZmZ2LVrV7Xvqs7vmUio48LlcmHOnDnIyMiAzWbDueeei2effbZavj179mDKlCkYMGAAYmNjoWkaNm7cqGxHbcdVfbdLxX/+8x8MGDAAMTExSElJwYQJE5Cbm1st3yOPPIKrrroKrVu3hqZpmDBhQsh1MUzEoDNnBLNnz9YB6CdPngxbmQcOHNB37txZ7+2JjY3Vx48fX6t66sJPP/2kp6Sk6Pfdd5/+0Ucf6Z9//rmemZmpR0dH65dffrnu9XprLOP111/XAei33nqr/sUXX+jLli3TExMT9WHDhtW474oVK3QA+po1a/StW7fqLpcrHIel67qu5+fn61u3btUfeeQRHYC+bdu2kMvYtm2bDkBfsWJFte+ys7P1rVu3hqGlkUeo4+LWW2/VbTabvnDhQv3LL7/UH3roIV3TNH3u3LlSvpUrV+rp6en6H//4R3306NE6AP3LL7/024a6jKv6bJeKjRs36haLRb/mmmv09evX66+//rreunVrvWfPnnpFRYWUNyYmRr/44ov1O+64Q7darY1y7TNMuOBJ1BlCfUyi6kJTmESVlJToJSUl1bYvWrRIB6B/9dVXAfd3u916enq6Pnz4cGn7G2+8oQPQ//WvfwXcv2oSdfjw4ZDbHixVdYR7EnUmE8q42L17t65pmj5v3jwp72233abb7XY9Ly/P2ObxeIz0e++9p5ys1HVc1Ve7AnHhhRfq3bt3l34IfP311zoA/fnnn5fy0voa69pnmHDBj/POMH777TfcfPPNSExMRGpqKiZNmoTCwkIpj67reP7553H++efDbrcjKSkJ/+///T8cOnRIyufvcd7p06cxefJkJCcnIy4uDqNGjcKhQ4egaRoyMzNDbo+maSgtLcWqVauMR09DhgwJV3cEJDY2FrGxsdW2X3TRRQCA7OzsgPt/8803OHHiBCZOnCht//Of/4y4uDh88MEHtW7bkCFD0LNnT2zduhUDBw6E3W5HRkYGVqxYAQD49NNPccEFFyAmJga9evXCunXral2XPzZu3IgLL7wQADBx4kTj3FSdY3+P8zIyMnDVVVdh7dq16NOnD+x2O7p164a1a9cCqHx82a1bN8TGxuKiiy7y+4hx+/btuPrqq5GcnIzo6Gj06dMH7777bliPrSZCGRcffvghdF2vNgYmTpyI8vJy6byYTMHdbsMxruqjXSqOHTuGbdu2YezYsbBYLMb2gQMHomvXrtXaW9f6GCaS4NF8hvGnP/0JXbt2xerVq/HQQw/hzTffxH333Sfluf322zFt2jRcccUV+PDDD/H8889jz549GDhwIH777Tdl2V6vF6NHj8abb76JBx98EB988AH69++PkSNH1ro9W7duhd1uxx//+Eds3boVW7duxfPPPx/wGD0eD9xud41/Xq83yF6T+eKLLwAAPXr0CJhv9+7dAIDevXtL26OionDuueca39eWnJwcTJw4Ebfeeis++ugj9OrVC5MmTcJjjz2GWbNmYebMmVi9ejXi4uJw7bXX4vjx43Wqj3LBBRcYE7ZHHnnEODe33nprwP1++OEHzJo1Cw8++CDWrFmDxMREXH/99Zg9ezZefvllzJs3D2+88QYKCwtx1VVXoby83Nj3yy+/xKBBg3D69GksW7YMH330Ec4//3zceOONWLlyZY1tboxxsXv3brRo0QJpaWlS3qoxUZsxEI5xVR/tClQXLdu3vnDWxTCRhqXmLExTYvLkyXjggQcAAFdccQUOHDiAV199Fa+88go0TcM333yDl156CU899RSmT59u7PeHP/wBXbt2xZIlS7BgwQK/Za9btw6bN2/GCy+8gDvuuAMAMGzYMFitVsyaNatW7bn44othMpnQokULXHzxxUEdY6dOnXD06NEa882ePdvv6lggfvzxRyxcuBDXXXed338KlLy8PABAcnJyte+Sk5Nx5MiRkOr2V/6///1v9O3bFwDQr18/tGzZEk8++SQOHDiAVq1aAQBatWqF888/H6tXr8Y999xTpzqrSEhIQM+ePQFU9new5yYvLw/ffPMNWrduLbXtpZdewoEDBxATEwOgcgXy2muvxX/+8x+MHj0aADBlyhT06NEDX3zxhbGiMWLECJw6dQp/+9vfMG7cuICrGI0xLvLy8vye/9jYWFitVmOMhEI4xlV9tCtQXVVt8yU5OTmsdTFMpMGTqDOMq6++Wvrcu3dvVFRUIDc3F6mpqVi7di00TcNf/vIXuN1uI19aWhrOO++8gKqcTZs2AQBuuOEGafvNN9+snETV1J7a8Mknn8DhcNSYr2qSESxHjhzBVVddhbZt2+Lll18Oej+VSq2u6rX09HRjAgVU/kNq2bIlMjIypGPr1q0bAAQ1gahvzj//fGMCBYi2DRkyxJhA0e1VbT5w4AB+/vlnLF68GACksfnHP/4Ra9euxf79+439/NFY4yLQea7LGAhmXNF+AgCz2Wx8H+52eTwe6LpufDaZTNKktr6uA4aJZHgSdYbRvHlz6bPNZgMA47HJb7/9Bl3XlROYjh07KsvOy8uDxWKp9osz0GSopvbUhu7du0s3cxWhxF4cPXoUQ4cOhcViweeff+73V7UvVceWl5dXrQ/y8/ODKiMQ/va3Wq3VtlutVgAI6bUM9YWqbTW1ueox8owZMzBjxgy/Zdf0CojGGBfNmzf3+wqI0tJSOJ3OWo2BYMfVkSNH0KFDB+n7L7/8EkOGDKmXdl1++eXGDykAGD9+PFauXCm115dwXAcME8nwJOosIyUlBZqm4auvvjImNBR/26po3rw53G53tRtjTk5OvbRVRbgf2xw9ehRDhgyBruvYuHEj2rRpE1Q7evXqBQD46aef0L17d2O72+3Gzz//jJtvvjmocpjKcQkAs2bNwvXXX+83zznnnBOwjMYYF7169cLbb7+NnJwcKf7op59+AgDjkWgoBDuuWrVqhW3btkn7VvVRfbTrxRdfRHFxsfG56pxVlfXTTz/hj3/8o7TPTz/9VKu6GKapwJOos4yrrroKTz75JI4dO1btsVxNDB48GAsXLsQ777yDO++809j+9ttv16lNNpstpJWpcD62ycrKwpAhQ+DxeLBx40a0b98+6Hb0798f6enpWLlyJW688UZj+/vvv4+SkhLlZKCpEI5Vw2A555xz0KVLF/zwww+YN29ercpojHFxzTXX4JFHHsGqVavw4IMPGttXrlwJu90eUHShIthxZbVa0a9fvwZrl2oS27p1a1x00UV4/fXXMWPGDJjNZgCVKsP9+/dj2rRpIdfFME0FnkSdZQwaNAh//etfMXHiRGzfvh2XXnopYmNjceLECWzevBm9evWSJkiUkSNHYtCgQbj//vtRVFSEvn37YuvWrXjttdcA1F663KtXL2zcuBGffPIJ0tPTER8fH3DVoeqXel3Jzc3F0KFDceLECbzyyivIzc2V3rDcpk0bY/Xh6NGj6NSpE8aPH49XXnkFQGX8ycKFCzF27FjcfvvtuPnmm/G///0PM2fOxLBhw2r1j6ohWLlyJSZOnIgVK1YEfFt0p06dYLfb8cYbb6Bbt26Ii4tDq1atQo4pCpYXX3wRV155JUaMGIEJEyagdevWyM/Px759+7Bz50689957AfdvjHHRo0cPTJ48GbNnz4bZbMaFF16I9evXY/ny5XjiiSekFduysjL861//AlA5wQAq4wxPnTqF2NhYXHnllQDCM67qo12BWLBgAYYNG4Y///nPmDJlCnJzc/HQQw+hZ8+e1V6zsGnTJpw8eRJAZZzV0aNH8f777wOo/KHWokWLGutjmIih0d5QxYQV1cstVS90fPXVV/X+/fvrsbGxut1u1zt16qSPGzdO3759u5Fn/Pjxevv27aX98vPz9YkTJ+rNmjXTY2Ji9GHDhunffPONDkB/+umna9WeXbt26YMGDdJjYmJ0APrgwYPr1BfB8uWXX+oAlH+zZ8828h4+fFgH4PfFgG+++abeu3dv3Wq16mlpafrUqVP14uLiGusP9LLNwYMH6z169Ki2vX379vqoUaOqbQeg33XXXco66Ms2n332WR2Avm7duhrb+NZbb+nnnnuuHhUVJfVJ1fmtbduq+nPRokXS9h9++EG/4YYb9JYtW+pRUVF6Wlqaftlll+nLli2rsa3hIpRxoeu67nQ69dmzZ+vt2rXTrVar3rVrV/2ZZ56pVm7VMfv7873OdL3246q+26Vi/fr1+sUXX6xHR0frycnJ+rhx4/TffvutWr7Bgwcr6wv1JZ8M09houh5EJCbDBODNN9/ELbfcgq+//hoDBw5s7OY0GapWhA4cOID27dtLLyqsK7quw+Px4LXXXsPkyZOxbds249HPDTfcgMOHD1eLp2EYhmFCgx/nMSHx1ltv4dixY+jVqxdMJhO++eYbLFq0CJdeeilPoGpJ586dAQAnT540gnXrykcffYTrrruu2nb99yDp119/PSz1MAzDnM3wShQTEmvXrkVmZiYOHDiA0tJSpKen49prr8UTTzyBhISExm5ekyIvLw+HDx82Pp9//vlhW406ffo0Dhw4YHzu3r279J4mhmEYpu7wJIphGIZhGKYWsHcewzAMwzDVeOGFF9C7d28kJCQgISEBAwYMwGeffWZ8r+s6MjMz0apVK9jtdgwZMgR79uxpxBY3PDyJYhiGYRimGm3atMGTTz6J7du3Y/v27bjssstwzTXXGBOlhQsXYsmSJXjuueewbds2pKWlYdiwYdJLWc90+HEewzAMwzBBkZycjEWLFmHSpElo1aoVpk2bZrzQ1eFwIDU1FQsWLMDtt9/eyC1tGFid54PX68Xx48cRHx/PxpkMwzBMQHRdR3FxMVq1alXrFw7XREVFBZxOZ1jK0nW92v82m80W0PILqHwx6nvvvYfS0lIMGDAAhw8fRk5ODoYPHy6VM3jwYGzZsoUnUWcrx48fR9u2bRu7GQzDMEwTIjs7O2jfzVCoqKhAc3scyuAJS3lxcXEoKSmRtgXyk/zpp58wYMAAVFRUIC4uDh988AG6d++OLVu2AKhuQJ+amhqUh+WZAk+ifIiPjwcA7Pvlf5WrUYp8piBWqYJdyKrLA1VvGJ/GBnNMoRKofR7yVTBPlb1BHKo3QJnemnev07kIZ/fVV7BiqKurpnAeUxBlNcTqbyRGMNBxG8w4p8fgO651RVmB9jHyKNun15gHQFD3S5rHTLaHOtYCDRXV9aOFWJ8qDy2nuLgYPc7pYvzvCDdOpxNl8GAcWsNaxzuDE168VnIM2dnZ0itpAq1CnXPOOdi1axdOnz6N1atXY/z48di0aZPxve8162+l60yGJ1E+VJ38+Ph4JCQk8CSqjjTEJEpX5OFJVHXqYxIVbIk8iVJDx62nMSdRirrDOokimSJxEkU3BzOJCrQtnFhhglWr453h95NXpbYLql6r1XghcL9+/bBt2zY8/fTTRhxUTk4O0tPTjfy5ubnVVqfOZFidxzAMwzARjlnTwvJXV3Rdh8PhQIcOHZCWloYNGzYY3zmdTmzatOmscq/glSiGYRiGiXBMGmCu4xzIBAReSvThb3/7G6688kq0bdsWxcXFePvtt7Fx40asW7cOmqZh2rRpmDdvHrp06YIuXbpg3rx5iImJwZgxY+rW0CYET6IYhmEYhqnGb7/9hrFjx+LEiRNITExE7969sW7dOgwbNgwAMHPmTJSXl2PKlCkoKChA//79sX79+nqLD4tE+D1RPhQVFSExMREHs08g3ueZcTDP0CmqWB0A8NAYA0UMgzpPcPEJNbUPkGMVVPFLdRkggeqm5cpxCP73Uq1E02fSwcQ9BSKYmBApv5RHpD0++T2KhtFsoa600/EY7L6moCOYApcr1S1tl3egH03Sdv8Fq2NQ/Jcj51EfGy1XdZ7cJBjJo7jego0bpPu4ybn3khFDxwRtR4XL63c7xUUa6Ds26T70WF2KQag6Vg/Z2aUISPSt2+utuQ+l9iny0xWXKLM441HkRFrM8kiItojPNvJdTJTZb9pqEWVJ5ZK01ew/D21fcVEROrdrhcLCwnrxD636n3S3uT1sdYyJcuhePOc5Wm9tPRvhlSiGYRiGiXDMYXicZ645CxMiHFjOMAzDMAxTC3glimEYhmEinHCo68whPs5naoYnUQzDMAwT4fDjvMiEJ1EKil0ewOmR5u00QNRJoiAryBclTrcowyle019U4ZLLJ9+Vkn3KyXaPIpjTSgIoVeko4uEU5XPl0e/oVza6v9l/kCYtS35RXs0B6r75pDYp7g51eZmobxxtMGWpgnKDq0+9g1kRMa0KpKZ9GyUFt4pzER1l8pvft0za5TSwOaiuDfEt8b47qOqQAsWDCiZXvayRlhN6ULtKVyMVpQhEV90DAKDMJa5jlyIfHY/0eqP9IV2r5CDovp2T7VLdSdHiXyU9Dgepm45Hp8f/mFcFtWcXCg+3r47mS99d2LqZkW4eE2WkixziHldYIdKdSNtjyXimgewJVhIk7ikz0m5rnFS3k1zwv5WJOo4UVIj06XLfwwEg3/vS48TbuxOjxb9IVcB5sauuchamKcOTKIZhGIaJcPhxXmTCkyiGYRiGiXA01F0JxlOo8MOTKIZhGIaJcHglKjJpMq84cLvdeOSRR9ChQwfY7XZ07NgRjz32GLxe8Txa13VkZmaiVatWsNvtGDJkCPbs2dOIrWYYhmEY5kylyUyiFixYgGXLluG5557Dvn37sHDhQixatAjPPvuskWfhwoVYsmQJnnvuOWzbtg1paWkYNmwYiouLG7HlDMMwDFM3qtR5df1jwkuTeZy3detWXHPNNRg1ahQAICMjA2+99Ra2b98OoHIVaunSpXj44Ydx/fXXAwBWrVqF1NRUvPnmm7j99ttDqu8fmw7DGhOHfUcLjG1uhQqjV+fmRjqvxGGk2yTFGOnt+09K+yQTVcrJ3FIjXVok9o9NECqR6FihdDETmY4lCKVefLTYFwDsVqHeibdZ/G6n+6uUZVSNVkIUN1QF1CxGrrtFvDgmasFgU1g20Lrp9iGtRTnwCiXU2qNCieOr1qIKKJW9BYXWTctSWlIoVHS+36lUjar8KpUaVQiR7qt2vmhLqOrJ6a1ZqReMDVGwFiQ0n8tLVWr+66bnhfYHtfegY4huB+T+kXvEv22SS9EfdHspUc9SBZ7vMagsVug4d7jFToVEvUvHRGqc1Ugn2sS1ZCFa9TKf+xJV8anO66lyUd/7P54w0ruPFRrpbunCFmTEOS2NdLxNVH5FpxSpXHqOTyuOKdkujqOg3E3Sohw6bn4l6TKXSBc7xL0ZABykz+m56ZIca6QvaiOOaS+579KxU0iVhCRNr/tTZeLYyksa5kd65SSoro/zmHDTZFaiLrnkEnz++ef45ZdfAAA//PADNm/ejD/+8Y8AgMOHDyMnJwfDhw839rHZbBg8eDC2bNmiLNfhcKCoqEj6YxiGYRiGqYkmsxL14IMPorCwEOeeey7MZjM8Hg/mzp2Lm2++GQCQk5MDAEhNTZX2S01NxdGjR5Xlzp8/H3PmzKm/hjMMwzBMHeGXbUYmTWYl6p133sHrr7+ON998Ezt37sSqVauwePFirFq1Ssrn6+Su63pAd/dZs2ahsLDQ+MvOzq6X9jMMwzBMbalS59X1jwkvTWYl6oEHHsBDDz2Em266CQDQq1cvHD16FPPnz8f48eORlpYGoHJFKj093dgvNze32uoUxWazwWazKb9nGIZhGIbxR5NZiSorK4PJJ1jXbDYbrzjo0KED0tLSsGHDBuN7p9OJTZs2YeDAgQ3aVoZhGIYJJ6YwKPNUFktM7WkyK1GjR4/G3Llz0a5dO/To0QPff/89lixZgkmTJgGofIw3bdo0zJs3D126dEGXLl0wb948xMTEYMyYMSHXt+dQPizRDkkJR9OUgzlCnVFWInyljhwW6pGT2bI6LydaKPc0MrLpo8dyovSjdZuIksRC1EleooqTyvS5cswW//VZiDqPqmlofVQx+MCwrkb69W1ZRjrRLhRFvsq3EqJ2KScKGqtKcUWOu4Kk/50ljpXWcV5avCjTJ4DgdIWoT/I4dBCfM6oII/1G1TtUYaVS2imGSjWoiEu1D1XkqfIEsrijaj2atip2Uqm7aN+UUx826pvo0xInUaA53G74Q3VjNynUjipFnu/5tpKCLYpgEjdRk6l8Ml3kmJoRL7VOSdFG+iei9PKFtvdUmbg/UBWribS1XaK4xtKJOu/waaE8pf3RMk5WwBaS8RxH6qaKt7vf/N5IV5QKpdkF3YUK747+bY10kVO0NYfcl3xVqFQVR9NUOUe98yrcIo/T7f+atlqCu5ioQjiP3IdPFFb4yy6pke1R/u991PN0+2HhE/hrTomRdleoz3044ZdtRiZNZhL17LPP4u9//zumTJmC3NxctGrVCrfffjseffRRI8/MmTNRXl6OKVOmoKCgAP3798f69esRHx8foGSGYRiGYZjQaTKTqPj4eCxduhRLly5V5tE0DZmZmcjMzGywdjEMwzBMfcPqvMikyUyiGIZhGOZshSdRkQlPohiGYRgmwuGYqMikyajzGIZhGIZhIgleiVLQpV0zWGPikJ4oFDgZKcKDKZ14wFElSVahMID61w/Ck4qq6ACggngvOYliTVbY+VfhRRHvquhYod5RqQepGq9yf3Ha2yQLleD/69PaSGeT47AplHPUV2p49zRxDNSXTmXeBV/FlUj3To33uz2aHAcV/VF1V6Dl7lhS1qEC0c/U3o+qvSS/O1IJVZzRLqeqnmoecv6t1KQ6qGJQBe1/2n/0Qvb1cXOTDfSXLP1RG6WQyNE8niC89nyJtYo2JkaLa4bu76voq4KqElXegLRNbp8+pqo6WVEp8tDj9koedGLfVnZZ/WaUT+ouIF50AEBEf5JKjUIVeXFWcQbpONpzssxIp8eJ/vtfvlCEHSuW1WdUAUvbRfvqmouE8o6q2sqJN+AzXwunB6okjCcKRadPpzsUCrv8UqHoo8o5J6nP6/GvLo4h/U/vV1af+11huSj3eL64f6UQRXFr4mdKj+PLo7lG+sQxYf1F77s2u8if2Ez8X3CV+z+/4caMMDzOC/K6ZYKHJ1EMwzAME+GYwvA4T2V4ztQefpzHMAzDMAxTC3glimEYhmEinLCo83ghKuzwJIphGIZhIpywqPP4cV7Y4cd5DMMwDMMwtYBXohRc1DEZ9th4XJCeYGxrGStUInaLQrnlEfnPTxPpdT8L9QcA/PJbMfxBVURUxULroKoS6itlJwof6gXVLEZWF7VIoN5cNvgjiShiaDBiIlH2dUq2G2mqcsoqFEqcNOL9BchqOwpV3tlUqjgPVZmJfakSiubxraqCqOqax5C+Uii3qO6IHh/1aLOQ7fnl1DdMVi3F28h5ogpASRUnzoWJqMM0r1BP5bv89188UcHRPgAAJ1FJ0e/cXv9p+muVqiubE3USVY1R5WG0j8+ZR/H8gO5DVYlU8WYj6iuqdqwgCjCXrxSRQNsSZ/WvaqTqSJfHf1l0ay7xmaOqu4vbJEr7UOGYw+2/3IMFQkFG/eRyS8V1f06KUJNlk+uKqjR7tRSqYUA+Vtr9dKjZIY4jxyHyb/210EjT41N5RB7Ok33jtu35zUj/sG6jkW7eqZeRpopiClURxyYI9ZudKKTpfa3cR/WYWyT6J/834W2XSpR0e46J4zvwY46Rpio8e7xon53cv1KSxP0uPVGkneJSqFf4cV5kwpMohmEYholw+HFeZMKP8xiGYRiGYWoBr0QxDMMwTIRj0rQ6v+eJ3xMVfngSxTAMwzARjmbWpDe516oMnkSFHZ5EKejfJhFx8QlIIEGaNChaFdj8S54IFk0jAYozLs2Qyi90iKBIOqxpIHU0qYQOfhrU/muxCBD9lQRW9m4pAlJP+dgSJBLbGFqWnRwftWygv15MCosOSkyUCLpMMMl2GPCKz5pOootJpLGukQDrolNG2lZyUuSxiWBa3UICT+l2m7CPAQBrNAmWJ4fhIk+16blQHalFYTtC728ZXllIoBGRgJe0S7eLgOSonH1iu0ucS8Q2M5Kp5cKSQrfFibRVnO9o3cf/hASme5LaiX1M4vIvcop96N7FDnF8pyvEOKJB2DQA2eWVIwRoEHfHJHGedBIo3r2F2E6D848WCjsTapxKYrDRvYU4bt8Y89RYcXw5Jf4tdWgAOA1k1yXrIpGnDbmmK8gXeeVy+T+fEgHXNCi+JQmqzikW55j2YRRpx9dZp430ha3FWDlN6tuVI4KoAaBbC3EN5BPbF9U5o8fXvYUYU9Qy5lSZGL/0fmAO8E89vXsfI+2oEGUV5uYb6dgkcUwJyaLuOCJ+iSfX7eadx4x07hERGA4Amlnc16wxog8+XfmBkXY7iZ1VXJKou805om4SiB5Dzlcn0jcDOiQb6bKSaLyE+sdk1iSroFqVwZOosMMxUQzDMAzDMLWAV6IYhmEYJtIxmyRT+lqhsQNxuOFJFMMwDMNEOJpJg1bHFz1p4Md54YYf5zEMwzAMw9QCXoliGIZhmAjHZNZgquNKlIlXosIOT6IUJNksiI+2SHYkVJVF1TvFRNk0ME7YuehEpaabZGuGclIWtSAxlwg1GjRyeoiSqtAtFEnpcUK50tlzQtRHVFxCg/J7sVRJRJRc1F5EwuPxn4fsS5VeyRZiM+P1UYrRZ/rkO80j1Du6RvKUFojspB2eXw+KfW1CTaORurW2PaSqLQ6hYnI3zxBlQbbFMbaT8AHJDoZst2niGKh601QoK6bcB38w0lFtOhlpqrCj6IlpRtqV1FbU9+sukSdKHLc3Rpxlr1Uea2XEdqS4grZXbC8lNjWnysQ5jidKTqrgo7ZAHYjq7rcSWY0pq7rE9kJSVmmRE/5ItIk6XGSsUJUfVbmmxcrnMbtItCUpWhwHFSjllVFrE9Ef1NrHSWxmyklfUvsearsDAIPaCtUZtc4pJ2W1IMovnapFFfYzFS567sTxdEqOkfO5fa6536GKvGNFQvnYklg/UQueQodIJxKF3AmiKnT61EXvi5R4YpkCap9C+qkgXyjnNr6y0kjHpWYYaVeZuK/FNG8l1dG6u1DYJaWK66pTz1FG+tQJcX9OSRcq2dF9RFm9Wort7RJF36juB8XF/vs73GimusdEaYrzw9QefpzHMAzDMAxTC3glimEYhmEiHH6cF5nwJIphGIZhIhzNzOq8SIQf5zEMwzAMw9QCXoliGIZhmAinciWqjoHlaJgg+LMJnkQpSLYBlfZNQs3gIkuh0bvXG2l7ufDJ0olSzBzfzEh7WwpFFgCkOcU+plOnxf5uom7yElWcVZSb7Kt4q9rXIdQtcAr1jWaxSvm0KKG08RSTukuF8qXiqFC/RXfqJtoaK5QrmokonshxS/jUTY8JZH/d7fS/3SGOQ6f9QesjfWZKbC7ylObJVROfOlN5oZGOIb5zHqKWo754ZUQZRdVdNnIFxeUK7zv3UZEGAM9J4fllSRP+dZ6W6UZa37XBSEe17SLSHtE3dHxQFV6JSRxDbpGskKOeiqlx4nzklop8nZKECqm1WYzNLLeo45e8MiNN1WEnitWKH5pv09HTor1E+UWVflQFFkd8Kz/7WXgR/qmX6DOqTo32kPEPoLNdnLNCTfRPiVOMo2ZEtZdg8/8PqkX5cSPtiWthpF1mUbdDoagDgPxy0Y5o8k/QSgwcnWSsURUwVQdTb798YiCYRBSUgOxnp2mijdSXkJ4XyfuQHEe7RKq6FGOQnpdDJ2UVqo1cEFQhV14ixmDfXkJ5elGGUJXO/se/jTS91qmvXftuov//MihDqvuyjsLPLt4q90kVTnJ8KuU1VdxqbtFuqramyuRYd8M80OGYqMiEJ1EMwzAME+FomgatjgbEmpcnUeGGY6IYhmEYhmFqAa9EMQzDMEyEYzKbYKpjTJRJ53WTcNOkevTYsWP4y1/+gubNmyMmJgbnn38+duzYYXyv6zoyMzPRqlUr2O12DBkyBHv27GnEFjMMwzBM3al6xUFd/5jw0mQmUQUFBRg0aBCioqLw2WefYe/evXjqqafQrFkzI8/ChQuxZMkSPPfcc9i2bRvS0tIwbNgwFBcXqwtmGIZhGIapBU3mcd6CBQvQtm1brFixwtiWkZFhpHVdx9KlS/Hwww/j+uuvBwCsWrUKqampePPNN3H77beHVqFmAjSTpM6IPXnASJfu32WkbR3ONdKmaNm3zCiOlAMA3mihFKO+ZxTzcaHw8paSiSDxT9LMRMlWIdRTnkKiTHP79yYDAF2h9ItKSBB5iNLPI6n+iE+dPdbvdnh9JrAK7zyqxlG2lagHabtNpG4vURtqRDUJACgR35lsRGkTLRSHSGgp8kQJRVc8UQNSpV6xW/yyi08TKkZLlFC7+baR4okW/ZzX+xojTRVW1JPPmiJUntuPi+M7VCDO96+nZZXan3oKNdSW7NNGmiqxDhaIW0FL4umWbBfnJZ14rFFFF411pX5wgKxGox57bqJGKyE+etSXjfq7FRM12hFyfFmFIo/LI4/lbi1En3cp+9lIJyYLdSRcYn9LTo6R1iuI6oyMTYuDqGpjmhlpq48K1WsX17QtTvQtPW7qw2cl4r7YKP+qsVziaZhK/DbjfdV5Zfnkg2hXQqJQy7WOF9cobVNBhTjW/xE1polIUl0kf5sk2bcvjqgrncnkvJJztp+MwYM54v7QLJ2oLrveaKQHnS987foTNV+X5vI1ZSZtpCo8cvnAbvG/EhNFFY3U05N6hVKVLOlXaA2zFhGWl23qoe0/f/58rFmzBj///DPsdjsGDhyIBQsW4JxzhGJS13XMmTMHy5cvR0FBAfr3749//vOf6NGjR4CSzxyazErUxx9/jH79+uHPf/4zWrZsiT59+uCll14yvj98+DBycnIwfPhwY5vNZsPgwYOxZcuWxmgywzAMw4SFqpiouv6FwqZNm3DXXXfhm2++wYYNG+B2uzF8+HCUloofE2f7E6AmsxJ16NAhvPDCC5g+fTr+9re/4bvvvsPUqVNhs9kwbtw45Pz+KzI1NVXaLzU1FUePHlWW63A44HCIX79FRUXKvAzDMAxztrBu3Trp84oVK9CyZUvs2LEDl156afifADVBmsxKlNfrxQUXXIB58+ahT58+uP3223HbbbfhhRdekPJpmrxcqet6tW2U+fPnIzEx0fhr27ZtvbSfYRiGYWpNOILKf38cWFRUJP3RhYRAFBZWvqQ4Obnyxab8BKgJTaLS09PRvXt3aVu3bt2QlZUFAEhLq4z7yCFxDQCQm5tbbXWKMmvWLBQWFhp/2dnZYW45wzAMw9QNk6bBZKrj3+8LCm3btpUWD+bPn19j/bquY/r06bjkkkvQs2dPAAj4BMj3f/GZSpN5nDdo0CDs379f2vbLL7+gffv2AIAOHTogLS0NGzZsQJ8+fQAATqcTmzZtwoIFC5Tl2mw22Gy2ats1r7vyjwQTeuJF0HH0iPFG2lRBAp7zT4hC4lNIHvn5sHvXF2IfGshIg7JjRdCxFExOglh1st1L7WfKSJuIdQogB4FL20ndpqQWfvN4SXC2ToNpy8TxUYsa6XggB6nLBZPAclou2S5Z4pBgeS/dl/SBb90ULUrUYW4uAq9NJCDcYxHHUU6CyWmaBn3rdMXTLAcam1IzxIey00Zy3YECI51DrDHyS8XxOd2ivrhoYqtB7FKo5UWzrsL6BgDe3P2bkaaB3v3bCHFDRqJobzSxI5ECa4llhMsk+qaE9Ifu436iWgSmwcw0wP2CdBH8/FOuOJfXEauXNGJdE2cVx2MJ8DbnXJsIhKVh2NRmxkOO1UQsXXSz6HN6HZscIq275FupiQSgm03iOysRMURbhbjBTX7P0i6kAd2tbWKce+i16pKvKdpGjXznaNFV1A1xrLQsEDelHi3FfcL3vFZxbop8L8kt9S9ioRZKDiIAoNu9pBJ6XtokiGvSRgKry1xyowod/sUpiSTwPs4q9qcx2rRuM7F0oUHj1CoKHnEvMjfB2J/s7GwkEPGQv/+Bvtx999348ccfsXnz5mrfhfoE6EyiyUyi7rvvPgwcOBDz5s3DDTfcgO+++w7Lly/H8uXLAVSexGnTpmHevHno0qULunTpgnnz5iEmJgZjxoxp5NYzDMMwTO3RzKa6GxB7K/dPSEiQJlE1cc899+Djjz/Gf//7X7Rp08bYTp8ApRN1ZU1PgM4kmswk6sILL8QHH3yAWbNm4bHHHkOHDh2wdOlS3HLLLUaemTNnory8HFOmTDGkluvXr0d8fHyAkhmGYRgmsgmLAXGI3nm6ruOee+7BBx98gI0bN6JDhw7S97V9AnQm0WQmUQBw1VVX4aqrrlJ+r2kaMjMzkZmZ2XCNYhiGYZh6JizviQpxEnXXXXfhzTffxEcffYT4+HgjzikxMRF2u52fAKGJTaIYhmEYhmkYqtTvQ4YMkbavWLECEyZMAMBPgHgSxTAMwzARTjhjooJFVykKaJln+RMgnkQp0NzOSqsW6bX/RJFRJlRVOlHyaMkiuM6bc1jkCWC9EkjNZuxPVGeeshz/28m7PrwuqqqSiTIRqxjSLnOiUHV5Ck6KNhHLE6qco3XTcqh60OXz8lLd499mxhQlhqI5WuxvSRfP4PUiYW2iu8S50GLE8Zhiya8fH0sbU4JQsCGFvA+MnmOFMieKLKOXl4sbS5GTWHHEkmOIby3V/ZtZ5DtUISwwHOQ8lRMrlRJic9KCqJPSiPXKgLZCXdfu1C5xCHHkOAEMzRDH2iyaKpWI8ogqlTSqBCV9S8Z/lC7a14wqFE2yBQlV4VErjjiyj91i8punewthKUItPegTjQCCPMk6hyoLaR2HS8X2jHihSNWJos4bK64LS0GW2E6Udr7qW80pLFM08o9IJ2NNMxGlHynLG0Usiahqj9xnzG5iv+SW1XleK1HMkbZHgagoyfGpRFRU1eYgClGqgky2y+e7QzNZlVoFtfah/W8jNiyxZExQ1Z7dS46V9Ks3UbbLohZM+eT6cUkKQJHfSgaShfSNRsa5qVgoW71xYnxYCsSrcMwlPvZS9YTJjDDERIWpMYxBk3lPFMMwDMMwTCTBK1EMwzAME+FoJg1aoKXXIMtgwgtPohiGYRgmwjGZQjcQrlaGhx8+hRvuUYZhGIZhmFrAK1EMwzAME+GE5T1RddyfqQ5PohSYC7JhdsfBG0f876h3ElFxUb88LUqhriNqMgAwxRNlicoXj3rCOYSii3rhOfJFm9wVQp1HVXAeH6We+7BoryVWKO+iYoXi0BpPfLOIwi6KWAV4yolShtThKhXllOeeluqmKrzYdKEcsnUQ3mZRbbsYaef/dhlpc5LwLtRdQg2Y/+23RrpZL2FSbe0xQKobpF0aOZe6jRxrlFCEaY4SURaxlmpG1IMny8RxFxHvrhYx8qXVwi4+J6QJfziqVDpWLM7Fr0XiXEaROIaUGFE32RWelp2M9AldtnNIt/pXGYZ6O9Wpao8omMhm6JDVWrQOquqqy72cKqxo2jfcg35H66OqrObkvORWiGsm1SbOEfXIK0kQSsfYslxRV7T8ThwTVeG5xLk0EVUvRSfbqVKPKvIkf0mqsI0i6lkAOmm7bhEDt5h4zUVRf0TSN9SbjqobbfQeRT0GS+XjoW2hCsBm0eK6KibXSTOPuA7NVH1Lxxq91xI/S61UlpolxIh7ajlRfFKFKIUOF5NL3F9B/FJp3ZKKl6pQNXnM1xdhecVBHfdnqsM9yjAMwzAMUwt4JYphGIZhIhzNZIJmquNKVB33Z6rDkyiGYRiGiXBM5jCo8/hxXtjhSRTDMAzDRDphiIkCT6LCDvcowzAMwzBMLeCVKBUVxYDFCzNVhkhKGaLUaEF82IqI51wboTizeGWFnDda+J5px38WaTNR3XiIGof668UK9ZWljVBl6aVCRectE4qi4l8OSnWXFQuvJ/rLxt5ctMlxWuwf2074wOlOoQykCkBb63ZGOpr4/8VmyJ5iKiydehtp50+bjbSJKPLM6eRYSR0tiOcfUjNE2i0rIhGfQr5zwB9UeURlS1r5aSPtsIhykogXHfX7yisn5w6Ah/intSbWaP8+Ls5Fu0RxjhNt4tJMjBbptsRHL8Emzp0DQpkkO+fJvnMmEKWS7t9IS6OegzSPrtiXjAPfX2UaUWipvOw0ojqTlFhU9kfSbtqM6s2vEdqKaOLdFkXiRU46qbecGGvNiXrtCMS4a2aRFVoJSUJhZ8k/ItpbSvwmy1Wea0L1Z04UY81rF9cnVUHSewkAaB7/Pp30WKkqlJ5XnXTOiRJxLbSLJn5y1DfUSn3+5O9oe80OcV7jiXoQJfT+Snz36Dgg5VO1oa8qkUJ9+Kh3HlVmRnlJP3n9e4167ERFTa4LT3yqSCO4e1xd0UxhUOdxTFTY4UkUwzAMw0Q4HFgemXCPMgzDMAzTJDh48CAeeeQR3HzzzcjNrVy1XbduHfbs2dMo7eFJFMMwDMNEOJUv2zTX8a9p/8vftGkTevXqhW+//RZr1qxBSUnlC5F//PFHzJ49u1Ha1LR7lGEYhmHOAqreWF7Xv6bMQw89hCeeeAIbNmyA1Spi6IYOHYqtW7c2Spuado8yDMMwDHNW8NNPP+G6666rtr1FixbIy8vzs0f9w4HlKkxmwGSGtzjf2KQlCKWMpCQh6AktRDpKKFe8vv5WVuHXZklOF3UQny2dqFu0whyRbiP84UwFv4o6qJ8WUaaVnTwt1Z3YSajtft2020iXHBPKwuY9OhhpZ57og+g2Qolo7SjabU4ix03q1qKJEgeyx1hpolD0WUrE8UV16y/2V3hXlSUL5aO5uVDtWSpOi33dskpJI/5YVNmjuYjiUPP/u6LYLlSCVqL2qSAqJ6rOc/r4dWXEiM+W/CwjfVlGZyNd5hLHmpEoxpdG1HU2TeTRdNEf1KesGrpQv6k8wiQUqjgo+kaClo/Q/flU+9Lji6LjnB63T/u8pASqjlRB1YNxxG+wlHjOFTlF/5+uEP1vt8hHqnlFP7gTxfWmJZBrhly7eoEY/1QB6y34TexL1MHeZqJMOgYrj4OMHYXPYAWRODqJujKZeAm2jhf1mYrJPyjSDl8vQK1E3CtMCs+7qFLikUfUdl5yTzT5jCOxg39Fqe93sRZxHF6yVmD2BOFjSMaUpNalY4gaDmp1GeXBYzKZYKpjYHhd929smjVrhhMnTqBDhw7S9u+//x6tW7dW7FW/NO0eZRiGYZizAH6cB4wZMwYPPvggcnJyoGkavF4vvv76a8yYMQPjxo1rlDbxShTDMAzDRDjhmAQ19UnU3LlzMWHCBLRu3Rq6rqN79+7weDwYM2YMHnnkkUZpE0+iGIZhGIaJaHRdx/Hjx/HSSy/h8ccfx86dO+H1etGnTx906dKl0drFkyiGYRiGiXA0LQwv2wwmrjFC0XUdXbp0wZ49e9ClSxd07NixsZsEgCdRaqLjAXsctNhmxiavTQRF09f+W06LANHShDZGOooEqlrzZOsVSMGVIuiSBqx7bSLQUm/di+QRgZ150SLg2f71fCMdN2iEkW43uYdUtadQ1E0vSmuCqM/aqr2RNhFbFTNJu1uIgO5iqwiCzysXwZg2n+Xjlk5haRH/m3g5mmP3FiNtSRUB5/q5lxjpMnOMkT5VSoN6RR2JNtEOqybbMXiJ3YTmLBNpEtxqcoh9PCQIuNAhgrNpGG8ysX2hAb5ti49IdeuaGDv0fEcXiHy2X/eLdqSRwEnFjVNlk+EbfCsH5/u3caE2OiCBtVIdJIA2YCC7qm6VhQwNFFfc5DUTCfBViDp8A8tNNECYXIser/8gcxM5Pi+JFY6z+g8cTomhwcsyOV4xVtMdQrChk2vanSLOsZn2/wlxr/AUnxZ5SD+ZrKJ8Z4xsvRJHrjkaNF5IrFdoIHxKFBErOAuNtCXviJH2EkspvUJcO4gW7QB8RCV2YU+lOUWguGTtYyaiDqtclrGdnFc67qqNQXr+SR3UukurENe35hFtpXVI5kROYs2Td0ykW4j7o6+Apb442x/nmUwmdOnSBXl5eY268uRL0+1RhmEYhmHOGhYuXIgHHngAu3fvrjlzA8ErUQzDMAwT4ZztK1EA8Je//AVlZWU477zzYLVaYbfLq7D5+fmKPesPnkQxDMMwTIRjMptgquMkqK77NzZLly5t7CZUgydRDMMwDMNEPOPHj2/sJlSDJ1EMwzAME+FoJq3u6jxTw7xdvb7IysoK+H27du0Cfl8fNNlJ1Pz58/G3v/0N9957r7HEp+s65syZg+XLl6OgoAD9+/fHP//5T/To0SNwYf7QdUDX4U7yf1LMJUJx44kXCrloCBWR5hQWA95YoWoDgLKoBPiDKrxSiMLL08y/JULy3k2iTQOuEAURZUxxSlepjtiYQ0a62WWkHVQlFU8sbogdQ0G02J5TIo51+6FTRjrOKsq5prmskPvFJBRvLVKEwjFueE8jXeoWappT5UIVd6pMKHkoCTYxjKl1SouYOH/ZAQAWomiKIvYpBTZhX5NArD+SzOK8xP+6w0h7k4QaMzohTVmfprCxoLYS7gKhXHT9esBIU0UkLMTSI0oounQXUUXZZIshjSi/TLFEMWUjaiiqnLMQNROxhpEUTCqVVCAJNa1DoaRS3eJ1lQbG6yEfPNJXmqJcTaEs9IKMf6IetEP0ralcqNcQI8YKteMBgFMVRPlIFHkgSq4CTWxPIkpQS5G4t1Do+QZRmSXjuHwcXqEE1ezNjHSqSbSp2C16xFJA/jEVijFINYxaFBl3NnHt0LFVWTk5B0T9pkSlPKWKPKKeBbXPMvvUTZDOPVHYSXYyXoXiU1KtivrMcc3ErgHUsPUFx0QBGRkZkg2WLx6PR/ldfdEkJ1Hbtm3D8uXL0bt3b2n7woULsWTJEqxcuRJdu3bFE088gWHDhmH//v2Ij49XlMYwDMMwTKTz/fffS59dLhe+//57LFmyBHPnzm2UNjW5SVRJSQluueUWvPTSS3jiiSeM7bquY+nSpXj44Ydx/fXXAwBWrVqF1NRUvPnmm7j99tsbq8kMwzAMUyd4JQo477zzqm3r168fWrVqhUWLFhn/+xuSJtejd911F0aNGoUrrrhC2n748GHk5ORg+PDhxjabzYbBgwdjy5YtvsUYOBwOFBUVSX8MwzAME0lUvbG8Tn9N+I3lgejatSu2bdvWKHU3qZWot99+Gzt37vTbWTk5OQCA1NRUaXtqaiqOHj2qLHP+/PmYM2dOeBvKMAzDMGFEM5thMptrzlhDGU0Z30UOXddx4sQJZGZmNtpbzIOaRF1wwQUhFappGj7++GO0bt26Vo3yR3Z2Nu69916sX78e0dHRyny+QWe6rgcMRJs1axamT59ufC4qKkLbtm3r3mCGYRiGYcJGs2bN/P6Pb9u2Ld5+++1GaVNQk6hdu3bh/vvvR1ycWu1Uha7rePLJJ+FwOGrMGwo7duxAbm4u+vbta2zzeDz473//i+eeew7791f6juXk5CA9XShdcnNzq61OUWw2G2w2W7XtjpZd4UhIoHZh0IhqwxMnlDmaWxyrJUsot9ztRVvhlbs6IXevka5IF8o0SfUh+ZYRtQpJnzxHeOTFRIntROSHZmWy2qcoUXh2xTT3b+JY6hYFlDhFm3ILxLF+dVS8HbZ9M6HYuc4u1D5HrbIy8p0dwn+qY4pQJ9mI/92f0oWqx0nUcgXl4mS0jBXKnNQY8evKXJRjpA+VyYpIK/UyJLEBZURE5PL6V3ckUy88orCi6jWq2KSKRgDQiarIGy38/aIOikfNutO/+lAJ9VKzk7Hio3gy2UVbNKt/n0Y9ilwDKhUdHYMKTz1fdZ6k6KOeaYrtqn0liKpKup36qO6kYml95HqlCi+N9Ce9UZfrIk++lmSkM47tNNLuZOGlBgCppHJ3jFDvWo+IVXRPqigLZPzrzcT9y0RUtl67GDeSkjZK/lEpKcoqxC/3Mou4fyd6iBderrheqbKTjhVan2pMAL7jgirsiJKOjjvS/5LKk6hQpe20HQFUcZKKUvKOdFfP7Asdz+Q+6k4R90rLSaGeNZUTf716hGOigC+//FL6bDKZ0KJFC3Tu3BkWS+M8WAu61gceeAAtW7asOSOAp556qtYNUnH55Zfjp59+krZNnDgR5557Lh588EF07NgRaWlp2LBhA/r06QMAcDqd2LRpExYsWBD29jAMwzBMQ8GTqMofNwMHDqw2YXK73fjvf/+LSy+9tMHbFNQk6vDhw2jRokXNGX9n7969aNWqVa0b5Y/4+Hj07NlT2hYbG4vmzZsb26dNm4Z58+ahS5cu6NKlC+bNm4eYmBiMGTMmrG1hGIZhGKZhGTp0KE6cOFFtQaewsBBDhw6N3PdEtW/fvuZMhMaKKZo5cybKy8sxZcoU42Wb69ev53dEMQzDME2aKoVdXctoyqhinPPy8hAbG+tnj/qnVg8RT58+je+++w65ubnweuXn0uPGjQtLw4Jh48aN0mdN05CZmYnMzMwGawPDMAzD1Ddn8+O8qvc/aZqGCRMmSHHMHo8HP/74IwYOHNgobQt5EvXJJ5/glltuQWlpKeLj46VZoaZpDTqJYhiGYRjmzCYxsVJUoes64uPjYbcLIZPVasXFF1+M2267rVHaFvIk6v7778ekSZOMeKMzFUtFISxWLzSn8FryxgnfOPPpX420OznDSLs69BeFUP8mnx8ArtRzRF2lwnfObBeKsvwYEVdm8YrJqod4fDmJmuy308SXq0JIzv6QJvv0HSzwr5z0EPVhSoxQzbi9Yvt5R/8t0kRxcyxdvOT0jSzhJ3djM6KSAfBw3I/iQ0cRBHjIKcbSz0Ra6CwSx1Hmoj56oj9KiV/e6QpxrIk+oks3KbfEJVQ60eTXmdWs+U1TxZOJqJa0sgKSh4yVaPkRMlXuaQe+E+lk4bdnbiFeCSIpxSSVFFGZEQ8zSTFol+v2Eu82SZEned75V9jp5iB98YJAlz4QpR595y9d2TYpVGAK7ztfNIW/mXTc5PqhR2elqj3S8ESbOC+u1sJ2KstHoJUeJ+qIImrAYy37iDroYVOFHfFgpG2VjicYXzrIijy7V6g/TWTcalS9KanoFMo5leoOCE55R8cU3a46R0GOO6q6pOPW5CAnh7ZdofrTiTLQE5MssrhE/xWnizdnF8c2zAuaNZNW95WoJmpAvGLFCgCV3nkzZsxotEd3/gh5EnXs2DFMnTr1jJ5AMQzDMEwkwTFRwOzZsxu7CdUIeRI1YsQIbN++HR07+n+/EMMwDMMwTH3w/vvv491330VWVhacTqf03c6dOxV71R9BTaI+/vhjIz1q1Cg88MAD2Lt3L3r16oWoKHk59+qrrw5vCxmGYRjmLEczmaXH/LUtoynzzDPP4OGHH8b48ePx0UcfYeLEiTh48CC2bduGu+66q1HaFNQk6tprr6227bHHHqu2TdO0RnlPA8MwDMOc0ZjM8hvba1tGE+b555/H8uXLcfPNN2PVqlWYOXMmOnbsiEcffRT5+fk1F1APBDWJ8n2NwdmAbo6CbrZCJ3YapV4xAD3x4t1ZVtI9Nk188GrUKkEevGaPCIJ0xIhgcvrEmiofS0jwdDQJeLYTuwiXTUTAFpLY8Q2/yoHkV9pPGOn/2cVj2V+LRL4Sp5gMV7hF3SdbDzPSvVqKwOZvD4pA1SjSvtnfyoHlf8sTdjfe80cb6W2H84x0VkGZkU4kAe6DM0SQZ9t4sd1CgiU7xYq2fp8vWzzsPyUCTAe2bQZ/kKajggSix5KAVG+M2FcK9iXXieZrL1F62kia4ondBwkQtrQU71dT2qqYTH7zKK0xqn1HRhitg+angbwkODiQzYbIo/t89r+P0uqFjOeQ7WCqVULyqdouWXyIPCYauE22x5HsXrNQLhQ5ZMuevHKxf6t40c8OYqdUTq4r2m12GqxNHXg0cR6jA/wvpOc7xinGvImIG2iwtU6tZaj9DylHDgwnAeNmn38htQgI91uH9AW9roKzcNE84jsa4K6sj15XRERiKc410oVvPW2k4yc9aqSduvxIqd4wmapZOtWqjCZMVlaW8SoDu92O4uJiAMDYsWNx8cUX47nnnmvwNoXco6+99ppfXzyn04nXXnstLI1iGIZhGIahpKWlIS+v8gd3+/bt8c033wCodFXRfX/ANRAhT6ImTpyIwsLCatuLi4sxceLEsDSKYRiGYRiBZjaH5a8pc9lll+GTTz4BAEyePBn33Xcfhg0bhhtvvBHXXXddo7QpZHWe6rXrv/76q/FCLIZhGIZhwgjHRGH58uVGeNEdd9yB5ORkbN68GaNHj8Ydd9zRKG0KeiWqT58+uOCCC6BpGi6//HJccMEFxt95552HP/zhD7jiiivqs60MwzAMwzQQ//3vfzF69Gi0atUKmqbhww8/lL7XdR2ZmZlo1aoV7HY7hgwZgj179tRLW9xuNx5//HGcOCFiem+44QY888wzmDp1KqxW/zGe9U3QK1FVCr1du3ZhxIgRiIsTb8K1Wq3IyMjAn/70p7A3kGEYhmHOekymMKxEhRbBU1paivPOOw8TJ070+/994cKFWLJkCVauXImuXbviiSeewLBhw7B//37Ex8f7KbH2WCwWLFq0COPHjw9ruXVF00OIxvJ4PPi///s/jBgxAunp6fXZrkajqKgIiYmJOJHzGxISEuAilidOj/90HPFv8JD8tGN9H4A6ST46rF3U2YGcmiQixzlVJpQnDtIOqlKjCj6zz6v+CyqE8q5FjJhH0+N448ffjPQ5KWLCHE3UU6fKhCqlc7J4g328TeShaiQA2H5cWCTEWUXdQ1oLVdCBYtHejkni1wVVyhzSmxnpDjbRDlOpUPk9c1D+jZAWJ+q4qI149EwVebFRou3UTqaVjdiDlJ820hpRcWkVxSIdQDkkqeqIGkpW5FHLE4Uli2SFEuDmqLLsUKnwaFupQov6ZKpuGwEUfCobFnl/3W8euq/U7kAKMEl557+90nmiado+VR1EAeaOTZG+OkwsmOxRoh2tYkV/Hjwtxg691pvbRZ4WJmElRNunEQsSzdcCxkPUYpLFDdnfLfaRLIqo8o5a0VBqY/+j6k/lOAhirASLarwEcy0QWyDdIq5VelsrKipC6/Q0FBYWIiFBttgKB1X/k3577ykkxNhr3iFQWWXlSP3z/bVqq6Zp+OCDD4wFFV3X0apVK0ybNg0PPvggAMDhcCA1NRULFizA7bffXqe2+uPaa6/FtddeiwkTJoS97NoSUkyU2WzGHXfcgX379tVXexiGYRiGqUeKimS/P5vNBpvNpsjtn8OHDyMnJwfDhwvfVJvNhsGDB2PLli31Mom68sorMWvWLOzevRt9+/at5qHXGC/7DjmwvFevXjh06BA6dOhQH+1hGIZhGMYXLQyB5b+/r7Bt27bS5tmzZyMzMzOkonJycgAAqamp0vbU1FQcPXq09m0MwJ133gkAWLJkSbXvGutl3yFPoubOnYsZM2bg8ccf9zsTrI/lTIZhGIY5qwmjOi87O1v6Xx3qKhTFV62vUvCHg0h88XfIk6iRI0cCqFw2ox1V1XFs+8IwDMMwkUtCQkKdFzzS0tIAVK5I0Rjp3NzcaqtT9UFFRQWioxVxew1IyJOoL7/8sj7awTAMwzCMAs1kglZH25a67k/p0KED0tLSsGHDBvTp0wdApXPJpk2bsGDBgrDVQ/F4PJg3bx6WLVuG3377Db/88gs6duyIv//978jIyMDkyZPrpd5AhDyJGjx4cH20I+Iwe10we10oJ355VOTm8IhlxfIykY4imaxE9uW7uumh4hPyHVXpmMlOVO1TWCFUNm0ShKLLpPmv2+KjzqNKv9NEqbf9uFCX9SfecglWkd9uEWXFku0WKhQjdaVZZZXaqM7CN67AITphxymhNjqYLxRJ3x4TefJKRR+0bSbampgh2ro1X/y6GtZZKAYB+dzQLqHbqTpSUlc6SsQHotiRFF1kqV3yIANkJZDCsy4YRZ6E6obok19X+ZnRlWSlaqlmRZ6cx3+TfOuQRqSvuszfvtIxKB4VBFSNEaWfm9hWSeo8/2rAYLBUnJY+d44V57hEE7+WzcW/kVzCC5IqYx3SzUGMec1N0qTPqnnDUYWdVb4GqjCVC9cJqjpT5Vf2h+/2YFSNSh9D1Tgn5z5A+UpPxSCuK80l7jlaBQm6JmPeGycUmFGkHVHehvLOa/iXbZaUlODAgQPG58OHD2PXrl1ITk5Gu3btMG3aNMybNw9dunRBly5dMG/ePMTExGDMmDF1a6eCuXPnYtWqVVi4cCFuu+02Y3uvXr3wj3/8o2lMogDg9OnTeOWVV7Bv3z5omobu3btj0qRJ/MZyhmEYhqkPGuE9Udu3b8fQoUONz9OnTwcAjB8/HitXrsTMmTNRXl6OKVOmoKCgAP3798f69evD/o6oKl577TUsX74cl19+ufSG8t69e+Pnn3+ulzprIuRJ1Pbt2zFixAjY7XZcdNFF0HUdS5Yswdy5c7F+/XpccMEF9dFOhmEYhmEakCFDhgQ09tU0DZmZmSEr+2rLsWPH0Llz52rbvV4vXK6aV7Trg5AnUffddx+uvvpqvPTSS7BYKnd3u9249dZbMW3aNPz3v/8NeyMZhmEY5mwmHAbCTd2AuEePHvjqq6/Qvn17aft7771nxGU1NLVaiaITKKDydewzZ85Ev379wto4hmEYhmHw++O8OgaGhzGwvDGYPXs2xo4di2PHjsHr9WLNmjXYv38/XnvtNaxdu7ZR2hRyjyYkJCArK6va9uzs7Hp7DsowDMMwzNnN6NGj8c477+Bf//oXNE3Do48+in379uGTTz7BsGHDGqVNIa9E3XjjjZg8eTIWL16MgQMHQtM0bN68GQ888ABuvvnm+mhj46B7Ad0Lq1l0EdUExRCPNfrImFjiSUKQEqesSKHedvQrqqord4svqMKueYxQ3xDrPJglEYp/pR4AxJjFTvHEvy7Z7n8SbFP48NmIJI96CVLlWylkJVpchVAFpThKRdouyurXUbzAlSqHKD/mi76h/nw9WwqfP9omADheLFRZrRNEuVbat9FUsSPym7P2izY1bycKJcohr40cq28AqFfx/jSVcoiiUDNpxLtNUiZZ5Mta8p2TviHHShVyVG0XhMJKmSfAPiovO6k/FEG0KhVWNQUZ/azof6p4U/n2ycpKce3R+jQylgEAFnFumhUfMdI5zc4x0p1jRJ5yciuW3r8HcU3qUeIYNLdQs0rqNah94KivpNcmvyTZKJf47kkejyr/uWBVjKp9ghk7QaoxJfWnol2qPKayApGHHjfxFaQKPsmH0FmmaHiYaQR1XiQyYsQIjBgxorGbYRDyJGrx4sXQNA3jxo2D2105kKKionDnnXfiySefDHsDGYZhGOZsRzOZodVxElTX/SOF7du3G28H6NatG/r27dtobQl5EmW1WvH0009j/vz5OHjwIHRdR+fOnRET4//9IgzDMAzDMHXl119/xc0334yvv/4azZo1A1D5yqWBAwfirbfequYJ2BDUOsosJiYGvXr1Qu/evXkCxTAMwzD1iWYSweW1/Qv20XuEMmnSJLhcLuzbtw/5+fnIz8/Hvn37oOt6o7xoE6jFSlRpaSmefPJJfP7558jNza1mCHjo0KGwNY5hGIZhGH6cBwBfffUVtmzZgnPOEbGF55xzDp599lkMGjSoUdoU8iTq1ltvxaZNmzB27Fikp6fXm1szwzAMwzBMFe3atfP7Uk23243WrVs3QotqMYn67LPP8Omnnzb4rG/+/PlYs2YNfv75Z9jtdgwcOBALFiyQZqS6rmPOnDlYvny58Qr6f/7zn+jRo0fI9XnMNnjMNkmxQ9Ve1G+NKu+8JH8Lm+jecpe8YuckMr44ovSj6rf8cpHHRVb8MhKFeuRYsVCJpMWK+tykfIePSs1kEb9G3MSny07UdpKHn8n/RDkaou6oKKECKnEKFZE9Sl4+LoHwtou1+HfgpsojSTVDfL1S40S6zEX7X5Rj8fnRdX6qUCRFk2M1e4iXmpt4qZGlb09rMYY0p1DpVPMtq8rjqw5SraJLPnVUmRacL57f7bVRy0nqKTpealZfKT3L4NMPUj663ew3TzBqK7khvuo8cRya5JHn/xxDJwo+k/82SeeFetz5/MKn/nzeGOEXmWIV+3hMRCEK/+g69RsU7aMKMhP1dQQACymNXEt6lLjeJAUaUcBKvoIKpWSg861CUohq/q+ZcD5qknwvqYqSqFBNxCPPS1V41KOQegxSRSNVRJpq5Z4WOo1g+xJpLFy4EPfccw/++c9/om/fvtA0Ddu3b8e9996LxYsXN0qbQu7RpKQkJCcn15wxzGzatAl33XUXvvnmG2zYsAFutxvDhw9HaamQFi9cuBBLlizBc889h23btiEtLQ3Dhg1DcXFxgJIZhmEYJsKpazxUOF7W2chMmDABu3btQv/+/REdHQ2bzYb+/ftj586dmDRpEpKTk42/hiLkKfTjjz+ORx99FKtWrWrQgPJ169ZJn1esWIGWLVtix44duPTSS6HrOpYuXYqHH34Y119/PQBg1apVSE1NxZtvvonbb7+9wdrKMAzDMOGEbV+ApUuXNnYTqhHyJOqpp57CwYMHkZqaioyMDOkxDgDs3LkzbI0LRGFh5TJr1Yzz8OHDyMnJwfDhw408NpsNgwcPxpYtW5STKIfDAYdDLGEXFRX5zccwDMMwTOMxfvz4xm5CNUKeRF177bX10IzQ0HUd06dPxyWXXIKePXsCAHJycgAAqampUt7U1FQcPXpUWdb8+fMxZ86c+msswzAMw9QVfmO5QW5urt+3A/Tu3bvB2xLyJGr27NlB5Xvrrbdw9dVXIzbWv71AXbj77rvx448/YvPmzdW+81UL6roeUEE4a9YsTJ8+3fhcVFSEtm3bwuXV4fLq0Elwarnbv7VJc7v/gXmixO13OyDbxlD7FBokTYO7aSA0DZ62KixZrC4RK1Zikh+70u6INpFjcooAVZOVnDdimaGTi9BLLSVIW+OsIo8UtA1AIwGtHk0EulpI3aB2JtRWhbTJVCHqa24XeU5XiLa6vD4B9eS4aR/QIFs5AFmUZSn41Ui7m7WCX8gFrVtUocI+KIJp6XGHbJkRSDGrqq+h3x+jOiYaBKyyu1HYx0jnLkC5UhYavE6bZFbcGr2KMn0tSIj1jhTYTKxDSsg+9JopIsKMBLKd/v9zkXBWq6+4QWWvQ8cFPW5q+UMDyMl1Qcuk9wDNx75Hp5ZB5DtvEP+8lU5AtRCAa7QsxbmnAf+SLRANuifXoUYC+HWbsJdSiUvCDk+isGPHDowfP954NxRF0zR4PAp7rXqk3mQFt99+O/r374+OHTuGtdx77rkHH3/8Mf773/+iTZs2xva0tDQAlStS6enpxvbc3Nxqq1MUm80Gm82/PxvDMAzDMJHBxIkT0bVrV7zyyitITU2NiFcs1dskyneWGI7y7rnnHnzwwQfYuHEjOnToIH3foUMHpKWlYcOGDejTpw8AwOl0YtOmTViwYEFY28IwDMMwDYlmMkGro7qurvs3NocPH8aaNWvQuXPnxm6KQQO94KLu3HXXXXjzzTfx0UcfIT4+3oiBSkxMhN1uh6ZpmDZtGubNm4cuXbqgS5cumDdvHmJiYjBmzJhGbj3DMAzD1AEtDI/ztKb9OO/yyy/HDz/8wJOo2vDCCy8AAIYMGSJtX7FiBSZMmAAAmDlzJsrLyzFlyhTjZZvr169HfHw8GIZhGIZpurz88ssYP348du/ejZ49e1Z7O8DVV1/d4G1qMpOoYB4PapqGzMxMZGZm1n+DGIZhGKah0LS6v9U9AmKI6sKWLVuwefNmfPbZZ9W+O+MCy5s6JQ4PNId8QqhVC1XXUUWdhUjA4m0iT16ZXFbLWGI7QsZ1nJVeJOL0UKVZCamvgzPLSOeUtTfSqW6hJInXyqS6dQdR+sU2F+koO2pC8/ofpFQZ5TUJdYvZJ3+UV6iTlEvTRNkmKYTIDSQmipwLs0jnE8OHKB+7GhqESC18aP9bqBqH7OtqIZaPq6nAjIJEspraTbJlCeJGppQq1UKRJ5Vbs7pPabESjDIwgPWKso0hxk9K/R9IgRfMPxwSIyLtrtzXfx9UO99E3ScdHcmXQH5Eu0kfUPsls1co50q8YoA1O7pF7BxLVGYAPMntRJqoyCi0yz2a/38DVOFmIo+BpLPo0+eaatgGYY6hGsG68huf/QMpNY0vRFmmklNG2mtPNNLFRM1M7+cem7gXUVW0y9xAwiTNFIZJVNOOiZo6dSrGjh2Lv//97wEFYw1JvfVo+/btqy21MQzDMAzD1Ia8vDzcd999ETOBAmoxicrOzsavv4p35nz33XeYNm0ali9fLuXbvXs32rZtW/cWMgzDMMxZjq6ZwvLXlLn++uvx5ZdfNnYzJEJ+nDdmzBj89a9/xdixY5GTk4Nhw4ahR48eeP3115GTk4NHH320PtrJMAzDMGcv/DgPXbt2xaxZs7B582b06tWr2tOuqVOnNnibQp5E7d69GxdddBEA4N1330XPnj3x9ddfY/369bjjjjt4EsUwDMMw4UbT6h4Y3sQDy19++WXExcVh06ZN2LRpk/SdpmlNYxLlcrmMN3z/5z//MSSF5557Lk6cOBHe1jEMwzAMw6DyZZuRRsiTqB49emDZsmUYNWoUNmzYgMcffxwAcPz4cTRv3ryGvZsOsVYT4qwmlDiF+sRMZvFWotqgqr28cqFGM5OVU6rUA4CcUqG6oUoZWkc7uyirlKjUoogyxPPjLiMd10+8xb0oqiUpU6oaVrN/9Y+ESt2iUOdRouA00pbCY9J3VG3niRfBgdKzeuobR326PKJcKymnyCU6sAXx0XP6eOfll4tjSo31L3rQaV+RNpm8xAOQ+oiZFR55wfrXqZRwkrlfiKZiAZbsJa+4IHzqlOUq2+2rSvRfd1CPFRTt04PIExBV3eZg1HwK1V0A5LFNJJxkO73uLfR8u8T1GWMVY9bRaZDI4jPOoz3l8Ieb5PMG0Xj6Whk61KjqVQs41kS+kIcwrTtY9abKM1CVh/r8kXtLjF2o88wVRUa6zCKUjh6qbgyvOYcak0lSk9a6jDMAp9OJw4cPo1OnTrBYGvclAyH36IIFC/Diiy9iyJAhuPnmm3HeeecBAD7++GPjMR/DMAzDMOGDA8uBsrIyTJ48GTExMejRoweysipf8TN16lQ8+eSTjdKmkHt0yJAhOHXqFE6dOoVXX33V2P7Xv/4Vy5YtC2vjGIZhGIZhAGDWrFn44YcfsHHjRkRHRxvbr7jiCrzzzjuN0qZarYPpuo4dO3bg4MGDGDNmDOLj42G1WhETE1PzzgzDMAzDhAar8/Dhhx/inXfewcUXXyy9PLl79+44ePBgo7Qp5EnU0aNHMXLkSGRlZcHhcGDYsGGIj4/HwoULUVFRwatRDMMwDBNueBKFkydPomXLltW2l5aWSpOqhiTkHr333nvRr18/FBQUwG4XNiHXXXcdPv/887A2jmEYhmEYBgAuvPBCfPrpp8bnqonTSy+9hAEDBjRKm0Jeidq8eTO+/vprWK2yKql9+/Y4duyYYq+mh8Ojo8KjI94m1DRUqWcn3nmOCqFYo953sg+ejN0tvqM+btSHr1wTfUxLsjmLxfb4ZkY62kR84/IOGWlPXAupbq1cKFGoZ5SkaAlC6aSpfMuIss+TkC7towxsVHnLUVUWUUZZSBbqQaYRtaHdJntaNSPng6rwJPWPV6H2IsekuUX/Kf3TfBU70ndBmGQGpTSruc98ob/VVEo9KPKo1HnBquJUfmb0vCqh9anKCdavUOXnF6JiMKDnH1Wm0eMjdXjJ2VD6xhE/S+pbSVVxvupbSdFK9rGQ7Q436U/SVruz0EifRLyRdpPrItkuyvE98ybN/zGZlGO+bmNKIogxAqLCk6DXMb0Nkv636+LeUq5HkfwNJM87i1eiLrvsMqxZswbz58/HyJEjsXfvXrjdbjz99NPYs2cPtm7dWu29UQ1FyD3q9Xr9OiX/+uuviI+P97MHwzAMwzB1Qde0MKjzmubLNjdu3Ain04mBAwfi66+/RllZGTp16oT169cjNTUVW7duRd++fRulbSGvRA0bNgxLly41vPI0TUNJSQlmz56NP/7xj2FvIMMwDMMwDAD06tULq1atauxmGIQ8ifrHP/6BoUOHonv37qioqMCYMWPwv//9DykpKXjrrbfqo40MwzAMc3ZzFj/OA4Di4mLptQb+SEhIaKDWCEKeRLVq1Qq7du3CW2+9hZ07d8Lr9WLy5Mm45ZZbpEBzhmEYhmHCxFnunde1a1fld7quQ9M0v6FG9U2t3hNlt9sxadIkTJo0KdztiRgcHh1Wtw4TCZ20k2jmQwUiQLGLvcJIe62xRpqGNlp8QjCjSGC6i5RLNqPURe1kRGndNRH86SnINdK5FSJ/y+YdRUE+vz6oxQG1cVEGigeDKr/PdslVxSnsKbwxSeQLRbA23ZcGjipsWEyuCroLdLMIBtVIOKBUlkcEj0r9ROuIEr+GNGeZyGMN8J60YPpHcazKYHzpfKmrVgW1ayA3HEX7lLYewQaTB3HcGhw15glmeyALkrARZB00mFz6t0WthKRxbiZpEpwt2Z+I4/aSIHHfwOZyEKUF6SoLGSSx+/4jvmjXy0ietgsrJupn0ixa1FfhFtsT4HONkUBsSaQhXa+KwH4Vwd6jVOIPVeA3teAh+1qcJUbabRVWL6VE9BNN7tlmUwNNTM7ylaj3338fycnJjd2MatRqEvV///d/ePHFF3Ho0CFs3boV7du3xz/+8Q907NgR11xzTbjbyDAMwzDMWcygQYP8viOqsQl5WvrCCy9g+vTpuPLKK1FQUGAsnyUlJWHp0qXhbh/DMAzDnPWwd15kEnKPPvvss3jppZfw8MMPS+7J/fr1w08//RTWxjEMwzAMg8pHcaY6/jXRSVT79u1hNptrztgIhPw47/Dhw+jTp0+17TabDaWlpWFpFMMwDMMwDFA574hUQp6WdujQAbt27aq2/bPPPkP37t3D0SaGYRiGYShVgeV1/WPCSsgrUQ888ADuuusuVFRUQNd1fPfdd3jrrbcwf/58vPzyy/XRxkYh1eJEQpRTaR3ROUFsz3cLRV60V+SxUFGUSV6KpIIRK7ETcJF8To9Qg7Qj3iYHy4SCJmPA/zPS7nKqehHlmJwBVghDVEApt6vsQQJctLIiz7/CRanuooodoqLTKoQljm6Lk3bRiPJOgtRhLs0z0p7Y5mJfN1EhESWPbpHtj/zlqYbCPkWVJyhLlgColIyq+vRg2qfY1xelaJAeX6C+qqmO2vxTUCm8gjluD1GAUcWmrwLWTOuw+M2nqfpAcS5U9jgWH3UY/UzvM9R6xZwoxrYjUVgzJVYUiWa4hPLUdOKIkfakn+u3HYCPojUYOx8VijGreQKMZapkDGLM67p/BSUdsxaH6A9ikAW4xHmxkX6qV85ydV6kEvIonzhxItxuN2bOnImysjKMGTMGrVu3xtNPP42bbrqpPtrIMAzDMAwTcYQ0LXW73Vi1ahVGjx6No0ePIjc3Fzk5OcjOzsbkyZPrq40MwzAMc3Zzlj/Oc7lcGDp0KH755ZfGbopESCtRFosFd955J/bt2wcASElJqZdGMQzDMAwjqDIgrmsZTZWoqCjs3r0bWoQdQ8hnpH///vj+++/roy0MwzAMwzB+GTduHF555ZXGboZEyDFRU6ZMwf33349ff/0Vffv2RWxsrPR97969w9Y4hmEYhmHAgeUAnE4nXn75ZWzYsAH9+vWrNv9YsmRJg7cp5EnUjTfeCACYOnWqsU3TtEY1AKwPNK+numpIMQCTiArFZbIZabdX95v2xemhih2xvdRNVGNEcdPCLuorJD5WzWwij6m8QBRk9lGQ1cHHjSpuJKVRsAqyIFRWynK9/r3eqF+eN1q4eJt81Hiam3i0ERWRHiU873SLOH8molSiddM80lqu5Eno0x9mxaWm6E+54XVcwpc8wupwfSrbpynz6FoU/ELPd+1bVKvHE5rKSy0Y5SlVigVQLtbq2qihPpo2UZVagHNK72HlUfFi/5QORprem8wqj7uYZiJPYY4on455yNeibhP/4LzR8SRTEPcZBVKfBdqXegsGcS4okgKQpHXVNezTB/XGWW5ADAC7d+/GBRdcAADVYqMCPeYL1XNP0zTs3LkT7du3rzFvrV62yTAMwzBMA8IrUfjyyy9rtd/p06exdOlSJCYm1phX13VMmTIl6AWhkCdRwczMGpvnn38eixYtwokTJ9CjRw8sXboUf/jDHxq7WQzDMAzTpDhT/p/edNNNQRsY33PPPUGXG/Ik6uOPP/a7XdM0REdHo3PnzujQoYPfPA3BO++8g2nTpuH555/HoEGD8OKLL+LKK6/E3r170a5du0ZrF8MwDMPUlnAYCIe6fyT8P73++uuxcuVKJCQk4Prrrw+Yd82aNX63e33DK2qguLi45ky/E/Ik6tprrzVioCg0LuqSSy7Bhx9+iKSkJEUp9ceSJUswefJk3HrrrQCApUuX4t///jdeeOEFzJ8/v8HbwzAMwzB1phEe50XC/9PExEQj3imYx3Eqjh07htatWwfM88Ybb+CWW24JqdyQJ1EbNmzAww8/jLlz5+Kiiy4CAHz33Xd45JFH8Pe//x2JiYm4/fbbMWPGjAaXIjqdTuzYsQMPPfSQtH348OHYsmWL330cDgccDhFwXFRU5DcfwzAMw5wJ+P6fs9lssNls0rba/D+tD1asWOE3HSrDhg3D119/rVzcefPNNzFx4sT6n0Tde++9WL58OQYOHGhsu/zyyxEdHY2//vWv2LNnD5YuXYpJkyaFWnSdOXXqFDweD1JTU6XtqampyMnJ8bvP/PnzMWfOnGrbNXcFNJdCWQRZqSGp+MxiIFZ4iOrFRzngISt5bpqP/FBIixX1FzpEkFsSEdsVkNi3Zm5yYVCVjI+HVVBeakFsV2ohQlTDAADIcSjLDaIsa/5RI+1O6ajMRxV53qho8QVR9WjlhUbaXJovmprURuQnfRnI14t67OlRdrEPUQxqoXq6Bdu3wfz6DEaZptgujy+5HMmTTFWWQvEWjBJLqbQLgFLRR/wmlcWSPBodK747qM4T9XeDf7WpVB31hXSHrgKj96mKf8400o7bnzTSiV5SB/UGlDwixb1IcwvVq+6QfeO06JpVl167+Cemucr91qcc83QY1MKzUdquPEdUEazw7aPNCMb7MQxUvmyzbuq6qv3btm0rbZ89ezYyMzOlbbX5f9oYFBQU4PXXX8crr7yCXbt2KfO1bNkSI0eOxBdffFHt1Qhvv/02JkyYgAULFoRcf8hrgwcPHkRCQkK17QkJCTh06BAAoEuXLjh16lTIjQkXvlLHqseM/pg1axYKCwuNv+zs7IZoIsMwDMMEja6H5w8AsrOzpf97s2bNUtYbyv/ThuQ///kPbr75ZrRq1QoLFy7E4MGDA+Zfu3YtPB4PrrnmGrhc4ofAu+++i3HjxmHevHm47777Qm5HyCtRffv2xQMPPIDXXnsNLVq0AACcPHkSM2fOxIUXXggA+N///oc2bdoEKqZeSElJgdlsrjZLzs3NrTabrsLfMibDMAzDnKkkJCT4XQyh1Ob/aX2TlZWFFStWYMWKFSgpKUFBQQHeffdd/OlPf6px37i4OHz22We49NJLcdNNN+H999/H+++/j7/85S94/PHHMWPGjFq1KeSVqFdeeQWHDx9GmzZt0LlzZ3Tp0gVt2rTBkSNH8PLLLwMASkpK8Pe//71WDaoLVqsVffv2xYYNG6TtGzZskB4/MgzDMExTwqvrYfkLlkj6f/ruu+9i+PDh6NatG3bv3o2nn34ax48fh8lkQrdu3YIup0WLFli/fj22b9+OK664An/5y18we/ZsPPjgg7VuW8grUeeccw727duHf//73/jll1+g6zrOPfdcDBs2DCZT5Zzs2muvrXWD6sr06dMxduxY9OvXDwMGDMDy5cuRlZWFO+64o9HaxDAMwzB1QUfd3u5fVUYoRMr/0zFjxmDmzJlYvXo14uPja97BDz/++KORXrRoEcaNG4frrrsOo0ePlr4L1bou5EkUUPmMdOTIkRgyZAhsNltEPB+t4sYbb0ReXh4ee+wxnDhxAj179sS//vWvJvGSUIZhGIaJFCLl/+mkSZPw/PPPY9OmTRg7dixuvPHGkF+hdP7550uvYtJ1He+++y7ee+8945VNtbGu03TfFz7VgNfrxdy5c7Fs2TL89ttv+OWXX9CxY0f8/e9/R0ZGBiZPnhxSAyKNoqIiJCYm4tSBH5EQHy8rj0zk6Sd9eRdRXnltYpbsIj3r9MjdTKed1LvKRYqNItXR/VtYhJrGQ1RmZncFaavaL83X78ogGKWeNwglFcXksz0IhVY13zl/ZQVTji+Ktishebz2ZmJfV5mfzDK+SirJ02//ZlFFl4uNtMMqxo7NVUp2FuofU5nwRPTGNvdfeTjVQiqFXLD+Z6G+16YOvmq1IaiXDwajKvS5pqTzr1DnBeVjFsztmV7r8PGKk9SjIpj2VEwrI52si7FGx5fKz9JUQV5E6HsdUc86el+MayHSRJ2qGlMmp2gTvQcrFay+n4NVrhqNUtwbghiPRcXFSOnUE4WFhTXGGdWGqv9JWcdz6lx+UVER2rVKq7e21ifl5eV499138eqrr+Lbb7/FiBEj8Omnn2LXrl3o2bNnjfsfPXq0xjxA6K4sId+lnnjiCaxcuRILFy6E1Soukl69ehkxUQzDMAzDhA9d18Py11Sx2+0YP348Nm3ahJ9++gndu3dHamoqBg0ahDFjxijfVg5UPspr27Yt2rdvX+MfAOzZswdud3A/RkOeRL322mtYvnw5brnlFpjN4ldH79698fPPP4daHMMwDMMwNeDVw/N3JtClSxfMnz8f2dnZeP3111FWVoabb75Zmb9Pnz7Iy8sLuvwBAwYgKysrqLwhx0QdO3YMnTt3rrbd6/VK715gGIZhGIapL0wmE0aPHo3Ro0cjNzdXmU/Xdfz9739HTEyMMg/F6Qz+xbYhT6J69OiBr776qtpzw/feew99+vQJtTiGYRiGYYLgDFlIqhdatmyp/O7SSy/F/v37gy5rwIABsNvtNWdELSZRs2fPxtixY3Hs2DF4vV6sWbMG+/fvx2uvvYa1a9eGWlzEIwdPi+7SrdH+skuBiB6vCBy1mOQgUvrRahYf8itEELhJF9sTo4nFRAUJICeB5bqZWDN4/ds3BGqvkiCCLpWB2rr8xJjaUCiDloO0dqixHJ82mYvFS+P0KPGCVRqgrTmFDQUN/DWVn66xTYFsSjQS9E+DycujRDA53cOcL5aSPcnCLV0ZTC5VFqDPQg24VZUrCQEUdh2++Uwh3m5CtSSqRT7pqlSNc2/N9izVyg3GUiRUVPv69qvJ/zXm2vaZkbYOvdVIe6JEgLHmIvcWcty6TdhkeCzk5cQ+whFqjaIK2vdS6xxyAtzkXhHtIMHuZHx57cSANlBguUKcErQNlT/qOh7rSDgex50pj/NCZePGjfVWdsgxUaNHj8Y777yDf/3rX9A0DY8++ij27duHTz75BMOGDauPNjIMwzAMw0QctXpP1IgRIzBixIhwt4VhGIZhGD+EQ13XlNV5kUr9v4iFYRiGYZg64Q3TX1Pn9OnTePnllzFr1izk5+cDAHbu3Iljx441SnuCWolKSkoK+q3kVQfFMAzDMAwTLn788UdcccUVSExMxJEjR3DbbbchOTkZH3zwAY4ePYrXXnutwdsU1CRq6dKlRjovLw9PPPEERowYgQEDBgAAtm7din//+9+NYjrMMAzDMGc6uh7cC+xrKqMpM336dEyYMAELFy6UPPSuvPJKjBkzplHaFLLty5/+9CcMHToUd999t7T9ueeew3/+8x98+OGH4Wxfg2PYvvzyPRICGB3qVKFC7A3yNaFiiSO+LSaflTwP6XZq6RJD9tFIngKHWIiNs4o8VjdVk/lXT1VT54XJHkFCZcnio9DRFQotqY0qNaGiXFWZ1XYPRjGl6g+F7YU7UdhnmBxiu6Ry8kFqL1Er6marn9yA+X9bjbS3/XniC4UNkbnklLS/J4FIf0NVyKmUfsFay9S3jUug8RtMfUFYGiltWwKpU4MZU3VQdWmBlGhkXOgWoSKW1HLZe4zkiU6XGelkogK2UushamsjKWnl8aRSqLrM4n5Z7hZ9SITJcJD7oN0ivrB7xbVUZ9sX1f2LELAOY7toa1FxMVK6nFfvti/7jhxHfB3LLy4qQreMVk3S9gUAEhMTsXPnTnTq1Anx8fH44Ycf0LFjRxw9ehTnnHMOKirU9936IuS72r///W+MHDmy2vYRI0bgP//5T1gaxTAMwzAMQ4mOjkZRUVG17fv370eLFi387FH/hDyJat68OT744INq2z/88EM0bx7EO2wYhmEYhgmJs907DwCuueYaPPbYY4Y7iqZpyMrKwkMPPYQ//elPjdKmkF9xMGfOHEyePBkbN240YqK++eYbrFu3jg2IGYZhGKYeCIe6rqmr8xYvXow//vGPaNmyJcrLyzF48GDk5ORgwIABmDt3bqO0KeRJ1IQJE9CtWzc888wzWLNmDXRdR/fu3fH111+jf//+9dFGhmEYhjmr0RGGwPKwtKTxSEhIwObNm/HFF19g586d8Hq9uOCCC3DFFVc0Wptq9bLN/v3744033gh3WxiGYRiGYQJy2WWX4bLLLqs5YwMQ1CSqqKgopEj+4uJiSX7YJNFMgGZSqsYkLzTiX9dcFwotOITSxRuTJBVvIkqsKDNV/Ig0VXs1p2o0jzhtkpKEKmgCKX+oCixUlVow3neEav5ZKv81leInCGWTplIO+dYtnT+xj8rjy9cXzCiG+HeZnKV+8+gW/0q7ygJI20k7QHzHaN3ejn3FdqU3IPFsjEtR1x0MQZxLLYBHYcgEoaQKyvPMNw+9ZoJQ6gWlyqoNQfSVUm2nWnoI1D66j4mMeXJdlf74nZFOI16OB4qEmi+O+IO2iIkz0mbqI+mLidzziIK5xCF8+Dyk6USMLKmOy10ik90iPEF1j8tIVxsTqnMs+Teq1KY1q3Ll7SRNjrk+8eo6vHVciqrr/o3N1KlT0blzZ0ydOlXa/txzz+HAgQPS65gaiqACy5OSkpCbmxt0oa1bt8ahQ4dq3SiGYRiGYQR6mP6aMqtXr8agQYOqbR84cCDef//9RmhRkCtRuq7j5ZdfRlxcXM2ZASNynmEYhmEYJhzk5eUhMTGx2vaEhAScOnXKzx71T1CTqHbt2uGll14KutC0tDRERUXVnJFhGIZhmBrx6pV/dS2jKdO5c2esW7eu2su+P/vsM3Ts2LFR2hTUJOrIkSP13AyGYRiGYZSEwfalqT/Pmz59Ou6++26cPHnSCCz//PPP8dRTTzVKPBRQS3UewzAMwzBMQzJp0iQ4HA7MnTsXjz/+OAAgIyMDL7zwAsaNG9cobeJJlALdHAXdHAUQzzsqyNCIf5rJKfzrqArIE088y3x+QlBVF/VSo558JWYRg2YnMhazo0S0w1lG2uz/dOrWWJ8NpC3U088rFDRKNZRCUadUPwVSRan2V3mSBdOmQKol2s8aUS0p9pE2U9EfVdQRvNFCwSopBqsV5r+9tB26WSij1Aohsj0q2n8enzok375QFXZB+BUG9JALolxZSeV/c1C+dgBoP9NrV1IZhlOFRwnmnNHzIh2gYvx7/Q/IQN6YWgVRoZJ7S8xVk420e6N4XU2X1p2M9M9tBhtpB/G7axcr7iea2yFVXaEJVaqJPDui/qAJRIVX5CRt1UQem4WONf/nq9qiSojjWepbIrDTleNLcZ8I1Y+ylnihw1vHpaS67h8J3Hnnnbjzzjtx8uRJ2O32oGO16wueRDEMwzBMhKOH4XFeE3/DgURjeeX5wpMohmEYhmGaBO+//z7effddZGVlwemUV/x37tzZ4O0J2YCYYRiGYZiGpUqdV9e/pswzzzyDiRMnomXLlvj+++9x0UUXoXnz5jh06BCuvPLKRmlTrSZRX331Ff7yl79gwIABOHbsGADg//7v/7B58+awNo5hGIZhGPE4r65/TZnnn38ey5cvx3PPPQer1YqZM2diw4YNmDp1KgoLCxulTSE/zlu9ejXGjh2LW265Bd9//z0cjsrAwuLiYsybNw//+te/wt7IRsFkAcwWOQxPCuQlAZsns8VuCckiu8IyBqgMXPeHTgK9Y3URsO702o20O0oEdtpoucQ2RLcKKxpzUY5Uhyexld+6QQI7lYHiikBoKeAzUICnZMsSYlA7hR53sAGlXkUwOa2blkvPEc1ub2akXWYRrGs2ieNxaLLtS4ybiAFIMK5uI0H/Hp8AYX8EE1wfAKVFjpQptAByPZDAQBHErRpfyqBxTRGErbLrCFQuIbh2UFumIMsJwmZGLkBxLSnySONf97Ed8ZDxTG2JXGLcmai91GXjjXSBLsZzt9JfjXSRva3IQ4LMk33ua9EQn70mcf3YSaB4CQkm95ClkbxysW+zaHFMJiKqiSLjkYoFgAAB3iobI+lDEGIW0mfyfSKAxVMY4cByICsrCwMHDgQA2O12FBdXCrzGjh2Liy++GM8991yDtynklagnnngCy5Ytw0svvSS9UHPgwIGN8jySYRiGYZgzn7S0NOTl5QEA2rdvj2+++QYAcPjwYeiNtMwW8iRq//79uPTSS6ttT0hIwOnTp8PRJoZhGIZhCPw4D7jsssvwySefAAAmT56M++67D8OGDcONN96I6667rlHaFPIkKj09HQcOHKi2ffPmzfX22vUjR45g8uTJ6NChA+x2Ozp16oTZs2dXi8zPysrC6NGjERsbi5SUFEydOrVaHoZhGIZpanh1PSx/TZnly5fj4YcfBgDccccdWLlyJbp164Y5c+bghRdeaJQ2hRwTdfvtt+Pee+/Fq6++Ck3TcPz4cWzduhUzZszAo48+Wh9txM8//wyv14sXX3wRnTt3xu7du3HbbbehtLQUixcvBgB4PB6MGjUKLVq0wObNm5GXl4fx48dD13U8++yz9dIuhmEYhmHqj+uvvx4rV65EQkICXn/9ddx4442wWCqnLjfccANuuOGGRm1fyJOomTNnorCwEEOHDkVFRQUuvfRS2Gw2zJgxo5opYLgYOXIkRo4caXzu2LEj9u/fjxdeeMGYRK1fvx579+5FdnY2WrWqDJx+6qmnMGHCBMydOxcJCQl+y2YYhmGYSMfjrfyraxlNjbVr16K0tBQJCQmYOHEiRo4ciZYtW9a8YwNRq5dtzp07Fw8//DD27t0Lr9eL7t27N/ir1wsLC5GcLJRwW7duRc+ePY0JFACMGDECDocDO3bswNChQ/2W43A4DIUhABQVFQGotEjQLdHQTETxRvajFivmZKLUIMor7cQvYnt6V6leT0KakS52eMgXIl3hEaoPt1PU0dpG8tM2EUWeHiXUfN6YJDkfVbhQlYmJeh/4t3eh/gg6UbUp1UWBVEqS/UyIFjBUlaWqIljFmqR0EvtQFZ7mEvY6XmKxUu70rxqLiVKrtXR7othO6tZt5Boi50JSWFEFElXaBVCphaqKo2Nbc1WITApbIem8+CqkQlWpqVCVEyWukcBjLTTLIMlqhKomdV9NmP/2BWxL1S5S3aJcHWY/uX3QyXH7lispV0l7y4tEuiTPSJotQiDULE68BZreN2I1l5EudpOx6ZXDJcwlJ420O7G1SHuF6s9iFm2i10mJS/RHFFG6UgWf2UxUe+Ygx5bSqopaevm//0h5yP1K10WfeaNE39Qn4Xgc1xQf55177rmYNWsWhg4dCl3X8e677yoXRRrDP6/WbyyPiYlBv379wtmWoDl48CCeffZZPPXUU8a2nJwcpKamSvmSkpJgtVqRk5PjW4TB/PnzMWfOnHprK8MwDMMwtWPZsmWYPn06Pv30U2iahkceeQSaVv3HjKZpkTuJuv7664MucM2aNUHnzczMrHECs23bNmmydvz4cYwcORJ//vOfceutt0p5/XWsrut+t1cxa9YsTJ8+3fhcVFSEtm3bKvMzDMMwTEPj1XV4zsKVqIEDBxqvMjCZTPjll1+a3uO8xETx+EHXdXzwwQdITEw0Jjc7duzA6dOnQ5psAcDdd9+Nm266KWCejIwMI338+HEMHToUAwYMwPLly6V8aWlp+Pbbb6VtBQUFcLlc1VaoKDabDTabTfk9wzAMwzQ2lbYtdZ1EhakxjYDb7ca4ceOk8JtIIKhJ1IoVK4z0gw8+iBtuuAHLli0znk97PB5MmTIl5ODtlJQUpKSkBJX32LFjGDp0KPr27YsVK1bAZJKfhw8YMABz587FiRMnkJ6eDqAy2Nxms6Fv374htYthGIZhmMjBYrFg9erVyMzMbOymSIQc9fnqq69ixowZUoCf2WzG9OnT8eqrr4a1cVUcP34cQ4YMQdu2bbF48WKcPHkSOTk5UqzT8OHD0b17d4wdOxbff/89Pv/8c8yYMQO33XYbK/MYhmGYJk2VOq+uf02Zyy+/HBs3bmzsZkiEHFjudruxb98+nHPOOdL2ffv2wRukj1WorF+/HgcOHMCBAwfQpk0b6buqV72bzWZ8+umnmDJlCgYNGgS73Y4xY8YYr0AIlUp1ng3wKOKpiArJS9RW2rGfSR4yR/XKijqnR6yrElEKylz++1B6Fk5ivHTi26RbxGPJAocop1m0PIk0eYiaRKVaUqjlZCUWVfORNum1UEnRqlXeVQrvO+Uh+HyhWxRKLuKR54gTj34d5BzFxBBFnuIcUUw+cXgO4ndIsTmLxT5lBaSx5Pii40Va5R1mUviqAdKdU/P9zg+amyiuFL5xqjFR7WkB7WeVD1+oPogqJadJVrVJXmpBKjUNyPWqWxT9p2q3b30hercF1z56HfooAyHuA/QbjR6TQ/hyaqfFj1ETaWvpujeMdNYfHzDSsVGk5R5ZmSaNT3JMzXUxznWzuJYKvOLaSyJ+efT+SJJwkedRVMEHACb4H6uq8akFsYagqR6fSR6gYVKg1sDZqs6jXHnllZg1axZ2796Nvn37IjZWvq9effXVDd6mkCdREydOxKRJk3DgwAFcfPHFAIBvvvkGTz75JCZOnBj2BgLAhAkTMGHChBrztWvXDmvXrq2XNjAMwzBMY+EJQ2B5XfdvbO68804AwJIlS6p9p2kaPB7/r/+pT0KeRC1evBhpaWn4xz/+gRMnTgCotIKZOXMm7r///rA3kGEYhmEYpr6edtWFkCdRJpMJM2fOxMyZM40XU3LMEcMwDMPUH17UXV0XeVOQpk+tX7YJ8OSJYRiGYRoCj1eX3t5e2zKaMo899ljA7+vLvzcQIU+iOnToEPDllYcOHapTgxiGYRiGYXz54IMPpM8ulwuHDx+GxWJBp06dmsYkatq0adJnl8uF77//HuvWrcMDDzzgf6cmiOZxQvM41WohogTyEvWbuYV427nmEi8F89pkFQGxj0Jz8q7PFLNQu1AlkNtGVv2Ij5vksUbal2Ql6im3/HIyTaXOU6lMiHqNolSKBVo01qiiz7+aRr0vUWVRHzf6nJxs1z0+SjRFHdRzMK/cf2AiVQVJToJEIWS1iG98VTAFFaJcqiryRgnlnZWMI6qkKiUK0VizQoFEKyOKTcBHwUb6RKnUk8aBf6UkbZ+vKk6Jaqyp1JjKNpHNAVSaUp+oFH2qMa9QEupSH1ClnjxulApTqbAgrj0V5gDHIF0nTrKZ3CtswluTKvWo52bcoBFGOjla9IeV3rx0+d7gJQphEN9FE/F/LG8m1MyxZOjQRRJaBd2uk+vKd1GFivVkn01xbjSQ86RaDFAFXyvOl+/9tb7Qw6DO05t4YPn3339fbVtRUREmTJiA6667rhFaVItJ1L333ut3+z//+U9s3769zg1iGIZhGEbGo8uve6htGWcaCQkJeOyxx3DVVVdh7NixDV5/2F5wceWVV2L16tXhKo5hGIZhGKZGTp8+jcLCwkapu06B5ZT3338fycnJ4SqOYRiGYZjf4ZdtAs8884z0Wdd1nDhxAv/3f/+HkSNHNkqbQp5E9enTRwos13UdOTk5OHnyJJ5//vmwNo5hGIZhGFbnAcA//vEP6bPJZEKLFi0wfvx4zJo1q1HaFPIk6pprrpEmUVUHMWTIEJx77rlhbRzDMAzDMJHP3Llz8emnn2LXrl2wWq04ffp0tTxZWVm466678MUXX0jWbFartXqBfjh8+HCYW113Qp5ERZqDcn2hOcugOc2ySkdSvRBVClXqxQvvNfOJfSKPj1JM+j1AyrXkZ4ldmrUSZXmIAoTUrfIwowq8gH5pCm86SX1Ct0t+ZP79/AKiUn6R8DzZ3yoILAqPNFOAkD+z/32o8qjURftTZI+LEuW6yS87L5EHOX0iOGlX0V+D5W5RR7GDePWROgodovI4G2kIHQdUcRlArSX1ZzCedSrVGBVomUOPCtA8Tv9fqJRsvkrLKkyKYwOU1y7dLl3HkoKPrraTMom6VLOQtI86Dyp1mOpaVG5X9EegsU2bQf0iiQrV5ON5V8XJVxYZ6eZ3zjbSp4lqtX0iGXdun3NPx5qzHP4w0x/hOvHzI6b2wRjlmn2883wVkkY7aN+qVJPBjHnp/irKNDmK/eUOO5H+OM/pdOLPf/4zBgwYgFdeeaXa9x6PB6NGjUKLFi2wefNm5OXlYfz48dB1Hc8++2yt6jx69ChKS0tx7rnnwhTkNRFuQq7VbDYjNze32va8vDyYzUFKnRmGYRiGCZoqdV5d/+qLOXPm4L777kOvXr38fr9+/Xrs3bsXr7/+Ovr06YMrrrgCTz31FF566SXD/UTFqlWrsHTpUmnbX//6V3Ts2BG9evVCz549kZ2dHa5DCYmQJ1Gq90w4HI6gl+QYhmEYhgmeqpWouv4Ble9Won8OR/2/62rr1q3o2bMnWrUST1hGjBgBh8OBHTt2BNx32bJlSEwU7xdbt24dVqxYgddeew3btm1Ds2bNMGfOnHpreyCCXoeviorXNA0vv/wy4uLijO88Hg/++9//ckwUwzAMw0Q4bdu2lT7Pnj273kN1cnJykJqaKm1LSkqC1WpFTk5OwH1/+eUX9OvXz/j80Ucf4eqrr8Ytt9wCAJg3bx4mTpwY/kYHQdCTqKqoeF3XsWzZMunRndVqRUZGBpYtWxb+FjIMwzDMWY7Xq8NbR3Vd1f7Z2dmS963NZvObPzMzs8YVnm3btkkTnED4s4zTdT2glRwAlJeXS+3dsmULJk2aZHzu2LFjjROx+iLoSVRVVPzQoUOxZs0aJCUl1VujIgHN64bmcUNXWIrIgar+Y8E8KR2MtDOmufRddNFxqa4qvNHCBkQn9h3UWoBaMyjtDdzCckFlYeGL5hbBvtTSRQpiJW3VFMGYeqD6dP/Bn1LQZogWMEHnCSJoOTlapG0k6JjGEtA0DW6N9ohA2qgoEcRbuY/YyW7x3/ZCj+ibYqdoqzmYoH2F2ADwtSdR2IUoRAnSuaTB1tL4CHC+FXYrSlsUVYA7FOPDLY7Nd9xp7jL4xUKvKyoWIdsV1zTNo4P2h09+ekgasaahgexkTEhnWBJHkHJocD2JvPY936qgc9o/epT4p0nrdhWJPttfLvKca84X+7rEUwjNJ+BfI0HWuiXaSLubtTHSZnJv0hWWUvS6ohEkJiI70XyD46UxpZhsqMYXuc9rXoXoIVRrnjDjDUNMU9UcLCEhQZqUqLj77rtx0003BcyTkZERVN1paWn49ttvpW0FBQVwuVzVVqh8ad++PXbs2IH27dvj1KlT2LNnDy655BLj+5ycHOlxX0MSsqzmyy+/rI92MAzDMAwTQaSkpCAlJSUsZQ0YMABz587FiRMnkJ6eDqAy2Nxms6Fv374B9x03bhzuuusu7NmzB1988QXOPfdcaZ8tW7agZ8+eYWlnqAQ1iZo+fToef/xxxMbGYvr06QHzLlmyJCwNYxiGYRimkkh/xUFWVhby8/ORlZUFj8eDXbt2AQA6d+6MuLg4DB8+HN27d8fYsWOxaNEi5OfnY8aMGbjttttqXBV78MEHUVZWhjVr1iAtLQ3vvfee9P3XX3+Nm2++ub4OLSBBTaK+//57uFyVS6c7d+6s8fklwzAMwzDhw6PrUlhAbcuoLx599FGsWrXK+NynTx8AlU+vhgwZArPZjE8//RRTpkzBoEGDpJdt1oTJZMLjjz+Oxx9/3O/3vpOqhiSoSRR9hLdx48b6agvDMAzDME2QlStXYuXKlQHztGvXDmvXrg1LfVOmTMFjjz0WtseNtSXkSLlJkyahuLj6G1pLS0ulaHmGYRiGYcJDlTqvrn9nCq+//nqNL+lsCEIOLF+1ahWefPJJxMfHS9vLy8vx2muv4dVXXw1b4xoTXTNVql5UijxF2kuVc0QNE+WRX2amExWel6qCHCVG2lRWIPITNY1kb6BSnFElTiAbEKLcohYRQSmBTP77QKk0CtRe1fZgFDGSmomcr0B2N+T4aB9o5YVGOk5lZ6KyEyHtMLtkZVgS6SuPSaiW6JPxpGih8PqtVLQpLZbUR9tNFZsWoraq1uf+b5ySWi4YywwFuh5AqReMglNx7lUKPqlIVX6f+lSWHdSyRrKiod2hsDrSVGpFn7ol1Z/Un16/eaS20uuTlu/1n8enWGUf6Ip7gJeUmxwtjslja2Gk6X3JF2+sUCFLxy2dNM3/dt2/7Y4UPaKywQF8jpXYyagsrGh+Mg50SY1cs1K42v21nvCg7uo8hTa6SaJ68XdDE/QkqqioCLquQ9d1FBcXIzpa/CPweDz417/+hZYtW9ZLIxmGYRiGOXvxeDzYvHkzevfuHVGvWAp6EtWsWTNomgZN09C1a9dq32ua1mivXWcYhmGYM5lIV+fVN2azGSNGjMC+ffuQlJTkN6yoMQh6EvXll19C13VcdtllWL16NZKTk43vrFYr2rdvL3niMAzDMAwTHiJdndcQ9OrVC4cOHUKHDh1qztxABD2JGjx4MIDKN5e3bdsWJsVbcRmGYRiGCS9erw5PmGxfmipz587FjBkz8Pjjj6Nv376IjY2Vvg/mLezhJuTA8vbt2wMAysrKkJWVBadTfkV+7969w9MyhmEYhmGY3xk5ciQA4Oqrr5aFHb/773k8DR86H/Ik6uTJk5g4cSI+++wzv983xkHUB3qUHbo1RumLR5UdHpPwf3ISJZvFLBRTlgDvJ6U+d0pPMarGCdEXr5p6R6UyVCn6PP4VWpLKiSj74FF4T/lClS9kZdPXj8uog6rlFPsqVX4BypX6gPoPUjWT1M8KhQ8tMoAq0eIgslyyv4mk23uFSlB3kl9btE3Em0zpS+ezj0qRJJVLxzxVihE1oKT4o/0aQCEnqd9UY42q2qRi/Pe/jgBKKoV6VJlWjB25ff4v5EB9rlJv0TElXaPygQdXH4UqLTWiJlQpJZ3iftL6vkeNdAGtWqXo9b3/KI7JS9Sj1C+PDiPatSq/vEAo/RzpuQBRsUp9rvLaUzSKlhlVs4I1HHjCsBJV1/0bm0i0nQt5EjVt2jQUFBTgm2++wdChQ/HBBx/gt99+wxNPPIGnnnqqPtrIMAzDMGc1PIkSYUWRRMiTqC+++AIfffQRLrzwQphMJrRv3x7Dhg1DQkIC5s+fj1GjRtVHOxmGYRiGYSIqnCjk6PDS0lLjfVDJyck4efIkgMqo+Z07d4a3dX5wOBw4//zzoWmaYXBYRVZWFkaPHo3Y2FikpKRg6tSp1TqZYRiGYZoaHq9Yjar9X2MfRd04efIkrrrqKsTHx6NHjx7o06eP9NcYhDyJOuecc7B//34AwPnnn48XX3wRx44dw7Jly5Cenh72Bvoyc+ZMv69S8Hg8GDVqFEpLS7F582a8/fbbWL16Ne6///56bxPDMAzD1Cd1n0DV/XFgY0PDiex2O9atW4dVq1ahS5cu+PjjjxulTbWKiTpx4gQAYPbs2RgxYgTeeOMNWK3WGs0H68pnn32G9evXY/Xq1dUC29evX4+9e/ciOzvbmGQ99dRTmDBhAubOndso0keGYRiGYcJDJIYThTyJuuWWW4x0nz59cOTIEfz8889o165dvbop//bbb7jtttvw4YcfIiYmptr3W7duRc+ePaVVqhEjRsDhcGDHjh0YOnRoSPW5LHa4LHaYiDpD9bZXp1uskZa7RR6dKDh8PY/irULNZrcLlZWb/FKgOiBaLiWayP5UAsAor/xIU3OVizYSrz9ZCRSEVx9Vt6g8qXwVcVQxpfCoUvrzhYqPsikodZ/KL0/lj6VSd/mopyTlEFEnSV5sVM0Uk+R3uwRV0UkKQ7VaSFLxSf5u1EMuCLWRwt+N+p/5lhuUD6IKhYpUwuxTt0qFJ5Wr1ZxHQTVlmvSlon+C8XFTDCmNeNEFUs+q7gPKc99MPEGgHow2i8jzW5m4B8TZhb1XrI/s2HL6VyNNx7CbHJRJobxzKR41yXpvUY7ZRzmtuj/rqgcuJtGfZoVXH125Ua3hOC0uxTfhhQPL/YcTde3atcHCifwR8iTKl5iYGFxwwQXhaIsSXdcxYcIE3HHHHejXrx+OHDlSLU9OTg5SU1OlbUlJSbBarcjJyVGW7XA44HAI6XYkuEIzDMMwDIVftinCiTIyMoxwooyMjAYLJ/JHUJOo6dOnB13gkiVLgs6bmZlZo9/etm3bsGXLFhQVFWHWrFkB82p+3uFS9RIuFfPnz2fPP4ZhGIaJcBoznEhFUJOo77//PqjCAk1W/HH33XfjpptuCpgnIyMDTzzxBL755hvYbDbpu379+uGWW27BqlWrkJaWhm+//Vb6vqCgAC6Xq9oKFWXWrFnSJLGoqAht27YN6TgYhmEYpj7x6GF4nNfEvfMaK5woEEFNourrLaEpKSlBHfgzzzyDJ554wvh8/PhxjBgxAu+88w769+8PABgwYADmzp2LEydOGMt669evh81mQ9++fZVl22y2apMzhmEYhokkOCZK4HQ6cfjwYXTq1Knew4lqos4xUQ1Bu3btpM9xcXEAgE6dOqFNmzYAgOHDh6N79+4YO3YsFi1ahPz8fMyYMQO33XZbrZR5FW4dUW4dNJzQRBbaLORDlNn/CpxHl719KA4Sae70+rfKoe/0UNXtIuVEk0BQM7WS8Qk01kmAKqWCDAfaXpcuAjhpO2xW8YHGl9L6aBC7L1JgriJAWxlUrSmsTFQWJ76fVTHEKjuNYALLAwQmSxYyJPiaWmBAYbeiKwLRg7ETAXwCjclYo+MgKj/LSLsT0vzuK40b0lZlwL4voQbn0zr0IM6x70p4HQLFdaW9i15jnsqMCrsoqSyFEEHyQiFjRTHWNF8BhbfmfqaiF5OjVOSJEiKXk2Xi2kuOJueCFulzfXviRdA5tXqhweRSLD8Jircqziu11aLB44FWVeiZoRMHX781o1zJ3YVsD2LO4fJVDdUTPImqfMnmPffcg1WrVgEAfvnlF3Ts2BFTp05Fq1at8NBDDzV4m+ogl4kszGYzPv30U0RHR2PQoEG44YYbcO2112Lx4sWN3TSGYRiGYerIrFmz8MMPP2Djxo2IjhYT/iuuuALvvPNOo7SpSaxE+ZKRkVFtZQeoXLFau3ZtI7SIYRiGYeoPt1eHuY4rSe4mvhL14Ycf4p133sHFF18srSp2794dBw8ebJQ2NclJFMMwDMOcTfDjvErbl6r3RFFKS0tDFraFizPmcR7DMAzDMGcuF154IT799FPjc9XE6aWXXsKAAQMapU28EsUwDMMwEQ6/bLPyvY4jR47E3r174Xa78fTTT2PPnj3YunUrNm3a1Cht4kmUAptFQ7RFU6ozqCCPZrGSLyQLlwBKIV1SnIjtUYpdqDpPV6hVvGahjLFYArzCgShioshyqJM0hB4rPT4zaQftJ6pLkmxlAMAThEUCVWWFrNpTqLggq8sktRxVOqnUebScoJR6gZRb/lVWtFwv0RdRawxJtUf6SVMouirzkcucKtDI/u6kNjW2STpuqhikKJSmlXWT86qyTKHnWFJr0Tz+rVqqlakqSwXNH8QivUaO1bd85ZhStUP3f46l+hRjNqAaU5WPjs+KYiPpjIo10i1IdvrP2ywpceW2us3CkqvU6f846P5m0iazqstJ35gCXFeS6o/uQwp2eoJok5meO/9qTHpvV6mzw41H1+v8nqem/p6ogQMH4uuvv8bixYvRqVMnrF+/HhdccAG2bt2KXr16NUqbeBLFMAzDMEyToFevXsYrDiIBnkQxDMMwTITDgeWRCU+iGIZhGCbCOZsnUSaTqUb1naZpcLsVYR71CE+iGIZhGIaJWD744APld1u2bMGzzz7r992RDQFPohiGYRgmwjmbV6Kuueaaatt+/vlnzJo1C5988gluueUWPP74443QMp5EKbGYNFhMmmxjJalS/CtzqALGQlQbXlkzIymuVGovOuCpNx0tS/aC8t9WBJqhK+q2EoWKSh0G4memXGj1UWupFEZSg2l7JT85cqx0Oz0XVLXnc2yax+m3iZpCLQdalsrnj2JS+6Wp1Ggq/zWTQhWkU+0jPd/KmtXQMSypKFX9T6DH46bZLbJqjx4eHc8m0mKqGKKqOLNZoT4k+F5XcuXU81Gt2gwJleLPp0zat0qlpuq6VHki0iyBjkGlqKRlkfPn2bvFSJdcKFSauaVi/KfHifNK1Wi6Jit/zR6HkbZbhBq2zCXa6yHnTKcqN+L9qbrXesn4l+5LgNxB9DyRMRxtUvU5vcGS+wlRE9OzGEXUqVGow3gKAY/uhSeQN2WQZTR1jh8/jtmzZ2PVqlUYMWIEdu3ahZ49ezZae3gSxTAMwzARztn+nqjCwkLMmzcPzz77LM4//3x8/vnn+MMf/tDYzeJJFMMwDMMwkcvChQuxYMECpKWl4a233vL7eK+x4EkUwzAMw0Q4Hq8O01kaE/XQQw/Bbrejc+fOWLVqlfI9UWvWrGnglvEkimEYhmEiHrcX0Oo4CXI30ZCocePGNZrBcE3wJIphGIZhmIhl5cqVjd0EJTyJUmByVcDkstacEbIiTNPLyRdUbeKjFFMpe0g+k8LryqxShwWrvFDsL3ne0bqVfl8hepP5lktVZ2S7siyVmol6spHNAdtEv3P7V+3B7P/yoOdb9vYL4Oul2MekUNtI3mjkzKhEjEovwQBIdRBfQV2h+JS8GZVdG0AxRfCSssykPq/iHEsKQI2qU1XtkBV9qv5XjzXF9akos9ru0jVWsw+fuqCax0dAFPlo2029hxppK/HDbJ8oxgT1zDQRxVo13z63UOdF6aKOxGD+01BFo+K4zQjQ5woFrNRe1bmnnp7B3NdoHnLM9cnZ/DgvkuFJFMMwDMNEODyJikzq8BOJYRiGYRjm7IVXohiGYRgmwuGVqMiEJ1EMwzAME+Gc7S/bjFT4cR7DMAzDMEwt4JWoGlAqcIJQ10kKjoCV1Ox1JSvTglCvBVTLefxvVilXglECqeoL5CenKNdLfKkkzzSFx53sl1ezJ2E16FUQzLtIVH5kQdYn+/6RYqnKU6rOfx9KSr1gFGeB2iidb1Gf00M8zxQ/Ys3UH8/XI5Jakumq7cRTT9H/ZrKD5EsX4HegNC6oHxr8ey0Go0gNRlXrm0/52z+Y8RLsGJb2CWIMU6VlbHMjHV+aa6Q98akiP+0/l1AgV6uJtpd6VQbh4afqT81RIrJEx6vLlO6j/v0m5dsJ9eEj0FNcX+eoFni8ep3fE8WP88IPT6IYhmEYJsLRdV0ybK5tGUx44UkUwzAMw0Q4Xq9e55gmjokKPxwTxTAMwzAMUwt4JYphGIZhIhxd1+v8OI4f54UfnkSp0L2A7lUHFgYZyKssXhEkHdzOiguhFjYsqmBO+bg9fvNIBGkOqepPGmytCn40aQr7ExqITgJVdZ82hXr/MKlCgmtxvtX9RvpcYX0jBVKrgmcVweq+ZXlJH9KAbieJGne7FMIDAg30lk6XTyernHComagUNE7399I+IEHiCgFFtTEfjBWOrqgj5GvSt9Nrvj8ELXwIFSl6WnkCSH4yPuzNRBZiZyIF4BOLoGrCG1eFSJtqPj5Vn9NrQRVM7tt/UlmK+wA9bC1oIU5koHvDEBPFj/PCDj/OYxiGYRiGqQW8EsUwDMMwEQ4HlkcmPIliGIZhmAjn9wiTOpfBhBd+nMcwDMMwTK05cuQIJk+ejA4dOsBut6NTp06YPXs2nE6nlC8rKwujR49GbGwsUlJSMHXq1Gp5mhpNahL16aefon///rDb7UhJScH1118vfX8mniCGYRiGqVLn1fWvPvj555/h9Xrx4osvYs+ePfjHP/6BZcuW4W9/+5uRx+PxYNSoUSgtLcXmzZvx9ttvY/Xq1bj//vvrpU0NRZN5nLd69WrcdtttmDdvHi677DLouo6ffvrJ+L7qBLVo0QKbN29GXl4exo8fD13X8eyzz4Zcn+b1VKo9grBkCUZl46tioZqZoIZ1XdRCvmu4QVi6KLV21M5C2peqDf0rf3zro4oyem2r7EV00lNaMJ3mc8NQWYqoVHhKyx+pjmDUjb7l1my741Go6Dz0tFJVlWSp4nucRNVIMrpJ2kXSKnGX6ogkQyJd/R21h6FtpCopWYUXxmcPHmKpoxzDpL4glGVBE+I1JlmhSF8E2R+S9VQQSj1F/wdjEVRNIRcV7TeftA9R94Eoa1XWRartge2N/B8rvfaCUXxK+yqOR/O4lO0IJ5EcEzVy5EiMHDnS+NyxY0fs378fL7zwAhYvXgwAWL9+Pfbu3Yvs7Gy0atUKAPDUU09hwoQJmDt3LhISEuqlbfVNk5hEud1u3HvvvVi0aBEmT55sbD/nnHOM9Jl6ghiGYRgmnBQVFUmfbTYbbDZbWOsoLCxEcnKy8Xnr1q3o2bOn8f8ZAEaMGAGHw4EdO3Zg6NChYa2/oWgSj/N27tyJY8eOwWQyoU+fPkhPT8eVV16JPXv2GHlqOkEqHA4HioqKpD+GYRiGiSSq3hNV1z8AaNu2LRITE42/+fPnh7WtBw8exLPPPos77rjD2JaTk4PU1FQpX1JSEqxWK3JycsJaf0PSJCZRhw4dAgBkZmbikUcewdq1a5GUlITBgwcjPz8fQO1P0Pz586XB1LZt2/o7EIZhGIapDeGYQP0+icrOzkZhYaHxN2vWLL9VZmZmQtO0gH/bt2+X9jl+/DhGjhyJP//5z7j11lul7zQ/j1l1Xfe7vanQqJOoYE+Q9/c3Fz/88MP405/+hL59+2LFihXQNA3vvfeeUV5tTtCsWbOkwZSdnR3+A2UYhmGYOuDV9bD8AUBCQoL0p3qUd/fdd2Pfvn0B/3r27GnkP378OIYOHYoBAwZg+fLlUllpaWnVFjQKCgrgcrmqLYA0JRo1Juruu+/GTTfdFDBPRkYGiouLAQDdu3c3tttsNnTs2BFZWVkAKk/Qt99+K+0bzAmqj2fBDMMwDNPUSUlJQUpKSlB5jx07hqFDhxqLHCYfccaAAQMwd+5cnDhxAunp6QAqY5ltNhv69u0b9rY3FI06iQr2BPXt2xc2mw379+/HJZdcAgBwuVw4cuQI2rdvD6AeTpDXXc17SwvV7y6QQo6W6w1CdWMKou5gt0u+UqEdE1W0SGoVhZLNV2lE+9Cr+1egURWXrO4SX1BFnUqJE0jNa9KDUOkodyaqQqI0clOBj88uHiI59Oj+fec8VJEURDNoP1mI/5zVLJ9HVV+ZqI8YbatKHUkaJasBRTrKZwhRXzzV+VMq8qiSSnUtBXmNhYrmUYznYFWyKnWfpJzzr+xUqkKDUNX65lMqaIPwDAxGIed7fUtDx1Lzj1Nd4XGnE0Wpzw4kHeieGszYCcKDUXWPo/dsT8O8RkfXw+CdV0+vODh+/DiGDBmCdu3aYfHixTh58qTxXVpaGgBg+PDh6N69O8aOHYtFixYhPz8fM2bMwG233dakhV9NQp2XkJCAO+64A7Nnz0bbtm3Rvn17LFq0CADw5z//GcCZe4IYhmEYJpINiNevX48DBw7gwIEDaNOmjVzn7xM3s9mMTz/9FFOmTMGgQYNgt9sxZswY4xUITZUmMYkCgEWLFsFisWDs2LEoLy9H//798cUXXyApKQnAmXuCGIZhGCaSmTBhAiZMmFBjvnbt2mHt2rX136AGpMlMoqKiorB48eKAk6Iz8QQxDMMwjNcrvyy3tmUw4aXJTKIYhmEY5mwlHLYt9RUTdTbTJN4TxTAMwzAME2nwSlRtUSlDJDWHQgEDBOXVJEE9psxE0WIhXlVBqMYAoJxscCqkWFaz5j9NvLFMbofYQeXR5aMiomocyROO7F7mEv1B67ZbiLrLUWykc7wxRjopWvRBfrnsUWcj+0eTcj266Fu6Wk7Vb1HkA1WyFTu8fvOXuuRzqlK8mYN4xxxVuElOb6SfvUQX5fGpm9ZBhXsmxfERZzMJHf5VhSpPPMBnKKh8y1RKqlC9C32eVQSjytJ0V835KVSNRq63auMcCtVrEP6KgTzh/BJQpaboQ6qKC+KYSr1iu5VcR9XuH/T0ByM6llZGyHnRqKqT5JfUpT6+fSqVp8qbkRJMn9P6pAupYf6N6t66C1DDaUfJVMKTKIZhGIaJcLxePQwxUfw4L9zwJIphGIZhIpxIfsXB2QzHRDEMwzAMw9QCXoliGIZhmAiHV6IiE55EKXBEN4MjOkFpb0HjJ2ngryp4VvfKQc7KIEgaeB1lF1lIwKfDTYKZSd1FJMi5hafUSOe6ZPuF5nZRR7OKXCPtiWtBcokDN1Wchj/0KBHQfdwh2kdjTaNMcqCxxSTamGIV6SiSx45y8cFFbGZKxXYafJsmugm6WxxbK1eB3F5d9INuSTTSpbp/iwl6HDE0uLWi0Egnk8B+nQT/N/OK/gcAU2me37Z7Ypsb6XyzaFM5CQ4/kC+OOzVWHEN+uQiKpoHhZp/g7pgocXw0uN5uEQvR1K6FBvOrzLvlYHX/4x/wsU8Jl41LIMEGRSXnDsbiQxUArrJb8Wmr5hZWILqJXt/EPkUVkByqpVSA/VVB4752LVW4ycOJ/ApxvqLNoi/dXrUwRVPYClFbIrpdEiiQL+h9QylO8IjxDyAoWxx6jQbTz7rUl+QeR/vAFoQ6JAx4dR1aHV9R4OVXHIQdfpzHMAzDMAxTC3glimEYhmEiHH6cF5nwJIphGIZhIhxdD8Mkih/nhR1+nMcwDMMwDFMLeCWKYRiGYSIc3avX+WWZ/Dgv/PAkSsHOE6WILTFJaiiqbKIKkxgibXITtYoXQs3h9sqLfi6v6HqqxIomiqk4UpbNIvLEkvo8ZHk2zkqUJBDtbqMTtRsArUQox0wVRUbaXCyUetQ2RrcK+ZsnsZWR3virKLfMJVQz56TEGunW8VR3B1iKckTdhadFmzxOhAJVF5ndFWI7UdNUs89wlYl8TtEH8QqlksqChCqvNJfoA91ZQfLLKjVYrH6/M7uFwijFJqxsaP+3TibqSo2UayMqM6dohychTao6t0KMkVSbSGeViv3JsIOXLFATFx3ZfoZIpjTJRiVI1VgwCrQgCPgvQWE/pFJc0ScdHsVjD6p8NHuJZYxLvsYkyxrNv0WIbiFqUZKHKuRUx+erwKSY6F5B+HyUeURZMeRSoF2QfORrsT2ti5E+bqGKXnknqsiT+9O/pQuF5jdTLbSkMPSvqvVpRlCKNJrDRe671CLLTc5piVP0a0mxj0qwnmAD4siEH+cxDMMwDMPUAl6JYhiGYZgIh9V5kQlPohiGYRgmwvF6dfntpLUtgwkrPIliGIZhmAhH93qqOV/UpgwmvHBMFMMwDMMwTC3glSgFA+OKkBCnS8oazSHUXZrLIdIKBQz1adLNPl1NFTtkfy/xT9NtiWQH8guCqMMsVNUWpOJJUttFx4vtRC3kMYm2FzlF3bmFQqFClYTtmwkPuWZE0mVyEcUaAI2o4iTlEPHB8pYSlRr1IHPKZRllWqP9phFl9ZP7d1xEYedVebT5/9XmJYo61S87zUc5pEXH+M/nFf5+WilJkzy6R1FHlELxR3z6ACCd+PN9c7qZkaaeeonRKpUaqU96EiA+0NZ5fW8p5EB0sg99qkCLlbYHoSQKpD2jVwP1dKPegA6Hx2+eBE2MD5dZXBdUkUf73BOdINVN/eXouaT9Sf3X4jwlRtpGrhlTGRkTvl5xvyP5wQHwkrZQpSb1byy0JovtkjeduL7bFu4XddiIjyepS6WuA2Q/RqrspPuovPYoVKlHL1Xf8aHqW1XdFKq2i7eJa6F5NNm38LhoK0kXlcg+mfUFr0RFJjyJYhiGYZgIR/d6wzCJqvmVF0xo8OM8hmEYhmGYWsArUQzDMAwT4egej/LRfihlMOGFJ1EMwzAME+HoehhionSeRIUbfpzHMAzDMAxTC3glSoU5CjBHSao6PUqodDQP8cbSJcmI//ICeF1J2ahHW8mpmncgdVOVIPWWg49STI+K9p+P7G8iv1jsRIXXhnjhtU8U6jCqeDK5iXKR9hMAb0ySyFchVHi07VoU8cJzKzz1SICk7FlHlI5lxXQPaJYokibKNhM1jlMoLckvQNWvQRNVBlpkxZRULtnfW0H8/GgdDnJMij7QiGKKqv9M8SlSvl1OocRqkyDOt9Usjrt5NFGhkvMnHYJZHJ/DTfo5gIiOKqaoYs2j2ElSWFH/NHL9pEVTYzTRZ5bT2XJZWftIPtLeUuEXaZPOK8lDy40WXpBafDMjbSLbzakdpLot1HuS3DdgVihG6bVrFefSQ9KSHx/NHyUrP6laj95PyqLFOCiuIMdH1GtOci2Ym3fzmyfGLPq/pUO+xkCvdzJevGYxVr2auB9RXzs6Iuj4INaiAcdQqvM3Iy3dW4iCmSqTqT9lHClHKyb3H6qqJvm9SW1E2uLTB/UEq/MiE55EMQzDMEyEw5OoyIQf5zEMwzAMw9QCXoliGIZhmAiHV6IiE55EMQzDMEyEwy/bjEyazOO8X375Bddccw1SUlKQkJCAQYMG4csvv5TyZGVlYfTo0YiNjUVKSgqmTp0Kp1MRmMwwDMMwTQSv1xOWPya8NJmVqFGjRqFr16744osvYLfbsXTpUlx11VU4ePAg0tLS4PF4MGrUKLRo0QKbN29GXl4exo8fD13X8eyzz4Zcn9cWL3lQAZBVeGav/+2q/L55QvxFQFU2ch0kqfDjg0euS/NSv7Ay8oXYn5ZlN/mfawelUPQ5bpXPoJeoZjSiLpKUZqQdsvqQKCip16Fve2ndCqWTnIduJ15odF+Fks1XESlBFWHUh4/4+UGh9NPMpFyiMPQ2a2Wkf41KlaprSdJURUkVV6Vu6vUmyqW+Y1HktEaTcqj6yddvUPqOjBel36TK/5Fud/m/bbmT2skbkjNqLNdEvRyJClLlU0dVtjU7+1WVRcYLTdPrUHF9S+NZcT/RneXSPp5EMRYOFBIVnlnUYSJXh5V0DR0TUYrt9Mh1CxmnADST/+OgbbdUFPrN44gWyl1N818f9b7zPUN5duETqEeLtMqPkfpFUvWnCf7vXyaX6Gd6Xryu4JTXzJlJk5hEnTp1CgcOHMCrr76K3r17AwCefPJJPP/889izZw/S0tKwfv167N27F9nZ2WjVqvIm8tRTT2HChAmYO3cuEhISAlXBMAzDMBELx0RFJk3icV7z5s3RrVs3vPbaaygtLYXb7caLL76I1NRU9O3bFwCwdetW9OzZ05hAAcCIESPgcDiwY8cOZdkOhwNFRUXSH8MwDMNEElWTqLr+MeGlSaxEaZqGDRs24JprrkF8fDxMJhNSU1Oxbt06NGvWDACQk5OD1FT5MUZSUhKsVitycnKUZc+fPx9z5sypz+YzDMMwDHMG0qgrUZmZmdA0LeDf9u3boes6pkyZgpYtW+Krr77Cd999h2uuuQZXXXUVTpw4YZSn+XkruK7rfrdXMWvWLBQWFhp/2dnZyrwMwzAM0yj8bkBclz+wAXHYadSVqLvvvhs33XRTwDwZGRn44osvsHbtWhQUFBixTc8//zw2bNiAVatW4aGHHkJaWhq+/fZbad+CggK4XK5qK1QUm80Gm81WbbtusUG32OSgYxL8KcXFKixWaPCz2ycK1enxH5ZqNvmf8NHNNAfNTwPGoQrQBdSB1MoAef9t1S2K/LRu3wlsoHb5a4fKRkeVn1bls1331KzU1FzUssbpNy21iETfKoP/fdtIbTlomgaTa6og+gAB67/TLFrO82uRCMGlgeWqIPMoxSmSxiAdd9L59jlftO30mPxXEdT4UPWH1yxb7dChQ+1F6KWnW4Thh4vqMhT5qdUIvfZsZnmcR5HPUbTfSHCyqdx/gDU9PmrhohMhgW4R9yyPWb5/VRBLnkSbOKgKciAOckPykrNB89D+Iw5BUhC2WZP/hXjJvxRpf3Jr8ugJJE0ylYnrp2WMKCcBwoZF08k1Vm2okPFFbZdUY4qWRQQe0nUczHXoM+7qC133VBNu1KoMJqw06iQqJSUFKSkpNeYrK6tUkJl8VGImkwne3xVUAwYMwNy5c3HixAmkp6cDANavXw+bzWbETTEMwzAMw4SLJhFYPmDAACQlJWH8+PH44Ycf8Msvv+CBBx7A4cOHMWrUKADA8OHD0b17d4wdOxbff/89Pv/8c8yYMQO33XYbK/MYhmGYJk3Vyzbr9scv2ww3TWISlZKSgnXr1qGkpASXXXYZ+vXrh82bN+Ojjz7CeeedBwAwm8349NNPER0djUGDBuGGG27Atddei8WLFzdy6xmGYRimbrA6LzJpEuo8AOjXrx/+/e9/B8zTrl07rF27toFaxDAMwzDM2UyTmUQxDMMwzNmK7vWG7HThtwwmrPAkSoGp/DRMUR61soMqaOh2MkapXsdi8ulqot5RqX8ouuItDXRfk2QRQbfLO2tUQQj/ai8tGFVcOKGKQekLasNCbDKCsXDx+KjlVDYbJC1vD6LuYCx/INu7SAobagFD5ceKm50WRZRARK1lKhcviY0vkt+Ldm6UsNTRbbHiCyc5bqpKJAoyEBUYtRfRE6iZDGm2NVb6fEJrZqSPFYs6XGTg2ixi3EYR8YisqKPjUfSNS5Jsy+pLr2IMS8o5Uh9tB1Ur2i0iHWc1k7TIb3aUyJU4iKJModoMZgzKimCSJvcTJ+Sx4vBQ5Z0giiowrSKt6ierdF5IU0keH0cpyUKG1kfFi5rmf7tJYb1ClYjwkvuVb7t11fXq8p8nmOuYoCnUvZrvua8ndG8Y1Hn8OC/s8CSKYRiGYSIcnkRFJk0isJxhGIZhGCbS4JUohmEYholwvF6P/ELlWsArUeGHJ1EMwzAME+HoHi+g1XES5RvExtQZfpzHMAzDMAxTC3glSoEnJgmeGPlN5+EVrInCJLs8hXce9fKi6haVp15tUCry6uK1F0D1oixLlfYG0Y4g65YbovDH0kh9RMWlB6Pa860iiiiMVCoi0CyKsuhyvNu/WsjrazJaXkq+zPW7j65IA8UiSW2XTh41kppZKKbMPt5+bcjntiaVhxnZR6WGpedI8urznz3gd/QAPYq0skyFL6TPscltJGmq6qXbiYJSeay0GW6h/osh3nKVn1Xt9VuUGrqv6t5XrUz/Mj5lf9A+DOK4dbPV7/ZA6MHeB/yhagdpt8fdMP9G2TsvMuGVKIZhGIaJcCL9jeVXX3012rVrh+joaKSnp2Ps2LE4fvy4lCcrKwujR49GbGwsUlJSMHXqVDidNRvDRzI8iWIYhmEYpk4MHToU7777Lvbv34/Vq1fj4MGD+H//7/8Z33s8HowaNQqlpaXYvHkz3n77baxevRr3339/I7a67vDjPIZhGIaJcHSvp+6B5fW4EnXfffcZ6fbt2+Ohhx7CtddeC5fLhaioKKxfvx579+5FdnY2WrVqBQB46qmnMGHCBMydOxcJCQmqoiMankQxDMMwTIQTzklUUVGRtN1ms8Fms/nbpVbk5+fjjTfewMCBAxH1u8vC1q1b0bNnT2MCBQAjRoyAw+HAjh07MHTo0LDV35DwJMqHqgDu4uJiP9+Frx5qtSAZWqhiu+sQWF7N9kWRseEDyxX5gkjXNbA81Lo1RYB7yHYw1dro/6YYVGC54lel7htYLu1fhyBbRWA4DSyHT2A5/aypAss1GliuGpw1Bx03CGENLKf7+w+uVx5rsH1Ql76qzb5BBGLXJbC8Vm2q58Dyqv8Ven1bZXlcyvj+UMoAgLZt20qbZ8+ejczMzLqWjgcffBDPPfccysrKcPHFF2Pt2rXGdzk5OUhNTZXyJyUlwWq1Iicnx7eoJgNPonyouiC6dO7cyC1hGIZhmgrFxcVITEwMe7lWqxVpaWnI2ftuWMpLS0vDDz/8gOjoaGObahUqMzMTc+bMCVjetm3b0K9fPwDAAw88gMmTJ+Po0aOYM2cOxo0bh7Vr1xo//DU/P5J0Xfe7vanAkygfWrVqhb1796J79+7Izs5uss9pVRQVFaFt27Zn5LEBZ/bx8bE1TfjYmibBHpuu6yguLpYeU4WT6OhoHD58OGwqNqvVKk2gAnH33XfjpptuCpgnIyPDSKekpCAlJQVdu3ZFt27d0LZtW3zzzTcYMGAA0tLS8O2330r7FhQUwOVyVVuhakrwJMoHk8mE1q1bAwASEhLOuBtDFWfysQFn9vHxsTVN+NiaJsEcW32sQFGio6ODnviEk6pJUW2oerzpcDgAAAMGDMDcuXNx4sQJpKenAwDWr18Pm82Gvn37hqfBjQBPohiGYRiGqTXfffcdvvvuO1xyySVISkrCoUOH8Oijj6JTp07/v717D4qy3v8A/l6UXZblIgLCIgiKyUWQFI+1eENTEEeD7HirBEajwQJEs8xjJXklVNLTxdRjqGWjTkFTx0Rx5BLeEsIRlUCPXKxZJBXBy+Gy7Of3h+NzeGQR3Mz12d/nNcMMz/N99nm+7/2o83Gf7+5Co9EAAMLCwuDv7485c+Zg3bp1uH79OhYvXoy4uDhJN+D8OVGMMcYYM5pSqURmZiaee+45+Pj4YO7cuQgICEB+fr6w3qpHjx7Yv38/rKysMHLkSMyYMQNRUVFYv369iWf/5/ArUQYoFAosX778kb7l80lhztkA887H2aSJs0mTOWd71AIDA3HkyJEuj+vXr5/oHXvmQEZ/+fsyGWOMMcbMD9/OY4wxxhgzAjdRjDHGGGNG4CaKMcYYY8wI3EQxxhhjjBmBmygDPvvsM/Tv3x9WVlYIDg7GTz/9ZOopPbSUlBTIZDLRj6urqzBOREhJSYGbmxuUSiVCQ0Nx7tw5E864cwUFBZg6dSrc3Nwgk8nw3Xffica7k6W5uRmJiYlwcnKCSqXC888/j99+++0xpjCsq2yxsbEd6vjss8+KjnkSs61duxZ/+9vfYGtriz59+iAqKgrl5eWiY6Rct+7kk2rtNm/ejCFDhggfMqnRaHDgwAFhXMp16yqbVGvGTIebqPvs3bsXycnJWLZsGUpKSjB69GhERESgpqbG1FN7aIMHD4ZWqxV+SktLhbG0tDSkp6fjk08+walTp+Dq6oqJEyca/OJlU7t9+zaCgoLwySefGBzvTpbk5GRkZWVhz549KCwsxK1btzBlyhS0PejLeh+DrrIBwKRJk0R1/PHHH0XjT2K2/Px8vPHGGzhx4gRycnKg0+kQFhaG27dvC8dIuW7dyQdIs3bu7u5ITU1FUVERioqKMH78eERGRgqNkpTr1lU2QJo1YyZETGTEiBEUHx8v2ufr60vvvPOOiWZknOXLl1NQUJDBMb1eT66urpSamirsa2pqInt7e/r8888f0wyNA4CysrKE7e5kuXHjBllaWtKePXuEY37//XeysLCg7Ozsxzb3rtyfjYgoJiaGIiMjO32MVLLV1dURAMrPzyci86obUcd8ROZTOyIiBwcH+te//mV2dSP6XzYi86oZezz4lah2WlpaUFxcjLCwMNH+sLAwHDt2zESzMt6FCxfg5uaG/v37Y9asWbh06RIAoLKyErW1taKcCoUCY8eOlVzO7mQpLi5Ga2ur6Bg3NzcEBARIIm9eXh769OmDQYMGIS4uDnV1dcKYVLI1NDQAAHr37g3A/Op2f757pF67trY27NmzB7dv34ZGozGrut2f7R6p14w9XvyJ5e1cvXoVbW1tHb5R2sXFBbW1tSaalXGeeeYZ7Nq1C4MGDcKVK1ewatUqhISE4Ny5c0IWQzmrq6tNMV2jdSdLbW0t5HI5HBwcOhzzpNc1IiIC06dPh6enJyorK/Hee+9h/PjxKC4uhkKhkEQ2IsKiRYswatQoBAQEADCvuhnKB0i7dqWlpdBoNGhqaoKNjQ2ysrLg7+8vNApSrltn2QBp14yZBjdRBshkMtE2EXXY96SLiIgQfg8MDIRGo4G3tzd27twpLJQ0h5z3GJNFCnlnzpwp/B4QEIDhw4fD09MT+/fvx7Rp0zp93JOULSEhAWfOnEFhYWGHMXOoW2f5pFw7Hx8fnD59Gjdu3MC3336LmJgY5OfnC+NSrltn2fz9/SVdM2YafDuvHScnJ/To0aPD/yjq6uo6/M9LalQqFQIDA3HhwgXhXXrmkLM7WVxdXdHS0oL6+vpOj5EKtVoNT09PXLhwAcCTny0xMRHff/89cnNz4e7uLuw3l7p1ls8QKdVOLpdj4MCBGD58ONauXYugoCBs2rTJLOrWWTZDpFQzZhrcRLUjl8sRHByMnJwc0f6cnByEhISYaFaPRnNzM8rKyqBWq9G/f3+4urqKcra0tCA/P19yObuTJTg4GJaWlqJjtFotzp49K7m8165dw+XLl6FWqwE8udmICAkJCcjMzMSRI0fQv39/0bjU69ZVPkOkUjtDiAjNzc2Sr5sh97IZIuWascfksS9lf8Lt2bOHLC0tafv27XT+/HlKTk4mlUpFVVVVpp7aQ3nzzTcpLy+PLl26RCdOnKApU6aQra2tkCM1NZXs7e0pMzOTSktLafbs2aRWq6mxsdHEM+/o5s2bVFJSQiUlJQSA0tPTqaSkhKqrq4moe1ni4+PJ3d2dDh8+TL/88guNHz+egoKCSKfTmSoWET04282bN+nNN9+kY8eOUWVlJeXm5pJGo6G+ffs+8dnmz59P9vb2lJeXR1qtVvi5c+eOcIyU69ZVPinXbunSpVRQUECVlZV05swZ+sc//kEWFhZ06NAhIpJ23R6UTco1Y6bDTZQBn376KXl6epJcLqdhw4aJ3rYsFTNnziS1Wk2Wlpbk5uZG06ZNo3Pnzgnjer2eli9fTq6urqRQKGjMmDFUWlpqwhl3Ljc3lwB0+ImJiSGi7mX573//SwkJCdS7d29SKpU0ZcoUqqmpMUEasQdlu3PnDoWFhZGzszNZWlpSv379KCYmpsO8n8RshjIBoIyMDOEYKdetq3xSrt3cuXOFf/+cnZ3pueeeExooImnX7UHZpFwzZjoyIqLH97oXY4wxxph54DVRjDHGGGNG4CaKMcYYY8wI3EQxxhhjjBmBmyjGGGOMMSNwE8UYY4wxZgRuohhjjDHGjMBNFGOMMcaYEbiJYmYnNDQUycnJZnXd2NhYREVF/alzeHl5QSaTQSaT4caNG50et2PHDvTq1etPXYt1LjY2VqjDd999Z+rpMMb+BG6iGHtEMjMzsXLlSmHby8sLGzduNN2EDFixYgW0Wi3s7e1NPRWzl5eXZ7Bh3bRpE7RarWkmxRh7pHqaegKMmYvevXubegpdsrW1haurq6mnAQBobW2FpaWlqafx2Nnb23MTy5iZ4FeimNmrr69HdHQ0HBwcYG1tjYiICFy4cEEYv3f76uDBg/Dz84ONjQ0mTZokerVAp9MhKSkJvXr1gqOjI5YsWYKYmBjRLbb2t/NCQ0NRXV2NhQsXCrduACAlJQVPP/20aH4bN26El5eXsN3W1oZFixYJ13r77bdx/7czERHS0tIwYMAAKJVKBAUF4ZtvvjHq+dmxYwf69esHa2trvPDCC7h27VqHY3744QcEBwfDysoKAwYMwAcffACdTieM//rrrxg1ahSsrKzg7++Pw4cPi25XVVVVQSaTYd++fQgNDYWVlRW++uorAEBGRgb8/PxgZWUFX19ffPbZZ6Jr//7775g5cyYcHBzg6OiIyMhIVFVVCeN5eXkYMWIEVCoVevXqhZEjR6K6urpb2bvKlZ6ejsDAQKhUKnh4eOD111/HrVu3hPHq6mpMnToVDg4OUKlUGDx4MH788UdUVVVh3LhxAAAHBwfIZDLExsZ2a06MMengJoqZvdjYWBQVFeH777/H8ePHQUSYPHkyWltbhWPu3LmD9evX48svv0RBQQFqamqwePFiYfzDDz/E7t27kZGRgaNHj6KxsfGB61kyMzPh7u4u3D57mNs3GzZswBdffIHt27ejsLAQ169fR1ZWluiYd999FxkZGdi8eTPOnTuHhQsX4pVXXkF+fn73nxgAJ0+exNy5c/H666/j9OnTGDduHFatWiU65uDBg3jllVeQlJSE8+fPY8uWLdixYwdWr14NANDr9YiKioK1tTVOnjyJrVu3YtmyZQavt2TJEiQlJaGsrAzh4eHYtm0bli1bhtWrV6OsrAxr1qzBe++9h507dwK4W5dx48bBxsYGBQUFKCwsFJrclpYW6HQ6REVFYezYsThz5gyOHz+O1157TWhaH6SrXABgYWGBf/7znzh79ix27tyJI0eO4O233xbG33jjDTQ3N6OgoAClpaX48MMPYWNjAw8PD3z77bcAgPLycmi1WmzatOmhasMYkwCTfv0xY3+BsWPH0oIFC4iIqKKiggDQ0aNHhfGrV6+SUqmkffv2ERFRRkYGAaCLFy8Kx3z66afk4uIibLu4uNC6deuEbZ1OR/369aPIyEiD1yUi8vT0pI8++kg0t+XLl1NQUJBo30cffUSenp7CtlqtptTUVGG7tbWV3N3dhWvdunWLrKys6NixY6LzzJs3j2bPnt3p82JoPrNnz6ZJkyaJ9s2cOZPs7e2F7dGjR9OaNWtEx3z55ZekVquJiOjAgQPUs2dP0mq1wnhOTg4BoKysLCIiqqysJAC0ceNG0Xk8PDzo66+/Fu1buXIlaTQaIiLavn07+fj4kF6vF8abm5tJqVTSwYMH6dq1awSA8vLyOs3dma5yGbJv3z5ydHQUtgMDAyklJcXgsbm5uQSA6uvrDY63f34YY9LEa6KYWSsrK0PPnj3xzDPPCPscHR3h4+ODsrIyYZ+1tTW8vb2FbbVajbq6OgBAQ0MDrly5ghEjRgjjPXr0QHBwMPR6/SOdb0NDA7RaLTQajbCvZ8+eGD58uHBL7/z582hqasLEiRNFj21pacHQoUMf6nplZWV44YUXRPs0Gg2ys7OF7eLiYpw6dUr0Ck1bWxuamppw584dlJeXw8PDQ7TWqv1z1d7w4cOF3//44w9cvnwZ8+bNQ1xcnLBfp9MJa4aKi4tx8eJF2Nrais7T1NSE//znPwgLC0NsbCzCw8MxceJETJgwATNmzIBare4ye1e5rK2tkZubizVr1uD8+fNobGyETqdDU1MTbt++DZVKhaSkJMyfPx+HDh3ChAkT8OKLL2LIkCFdXpsxZh64iWJmje5bS9R+f/tbPvcvcJbJZB0ee/8tos7O/SAWFhYdHtf+tmJ33Gvc9u/fj759+4rGFArFQ52rOxn0ej0++OADTJs2rcOYlZVVh+fyQVQqlei8ALBt2zZRkwvcbVLvHRMcHIzdu3d3OJezszOAu2uqkpKSkJ2djb179+Ldd99FTk4Onn322T+Vq7q6GpMnT0Z8fDxWrlyJ3r17o7CwEPPmzRNq9uqrryI8PBz79+/HoUOHsHbtWmzYsAGJiYndej4YY9LGTRQza/7+/tDpdDh58iRCQkIAANeuXUNFRQX8/Py6dQ57e3u4uLjg559/xujRowHcfcWipKSkwyLx9uRyOdra2kT7nJ2dUVtbK2o8Tp8+LbqWWq3GiRMnMGbMGAB3X5kpLi7GsGHDhEwKhQI1NTUYO3ZstzJ0xt/fHydOnBDtu3972LBhKC8vx8CBAw2ew9fXFzU1Nbhy5QpcXFwAAKdOnery2i4uLujbty8uXbqEl19+2eAxw4YNw969e9GnTx/Y2dl1eq6hQ4di6NChWLp0KTQaDb7++usum6iuchUVFUGn02HDhg2wsLi7fHTfvn0djvPw8EB8fDzi4+OxdOlSbNu2DYmJiZDL5QDQ4c8AY8x8cBPFzNpTTz2FyMhIxMXFYcuWLbC1tcU777yDvn37IjIystvnSUxMxNq1azFw4ED4+vri448/Rn19/QNfgfHy8kJBQQFmzZoFhUIBJycnhIaG4o8//kBaWhr+/ve/Izs7GwcOHBA1CAsWLEBqaiqeeuop+Pn5IT09XfRZQ7a2tli8eDEWLlwIvV6PUaNGobGxEceOHYONjQ1iYmK6nSspKQkhISFIS0tDVFQUDh06JLqVBwDvv/8+pkyZAg8PD0yfPh0WFhY4c+YMSktLsWrVKkycOBHe3t6IiYlBWloabt68KSws7+oVqpSUFCQlJcHOzg4RERFobm5GUVER6uvrsWjRIrz88stYt24dIiMjsWLFCri7u6OmpgaZmZl466230Nraiq1bt+L555+Hm5sbysvLUVFRgejo6C6zd5XL29sbOp0OH3/8MaZOnYqjR4/i888/F50jOTkZERERGDRoEOrr63HkyBGhOff09IRMJsO///1vTJ48GUqlEjY2Nt2uDWNMAky2Gouxv8j9C7yvX79Oc+bMIXt7e1IqlRQeHk4VFRXCeEZGhmghNRFRVlYWtf/r0draSgkJCWRnZ0cODg60ZMkSmj59Os2aNavT6x4/fpyGDBlCCoVCdK7NmzeTh4cHqVQqio6OptWrV4sWlre2ttKCBQvIzs6OevXqRYsWLaLo6GjRIna9Xk+bNm0iHx8fsrS0JGdnZwoPD6f8/PxOnxdDC8uJ7i7ednd3J6VSSVOnTqX169d3eD6ys7MpJCSElEol2dnZ0YgRI2jr1q3CeFlZGY0cOZLkcjn5+vrSDz/8QAAoOzubiP63sLykpKTD9Xfv3k1PP/00yeVycnBwoDFjxlBmZqYwrtVqKTo6mpycnEihUNCAAQMoLi6OGhoaqLa2lqKiokitVpNcLidPT096//33qa2trdPn4WFypaenk1qtFv7c7Nq1S7RYPCEhgby9vUmhUJCzszPNmTOHrl69Kjx+xYoV5OrqSjKZjGJiYkTXBi8sZ0zyZERGLOxg7P85vV4PPz8/zJgxQ/Qp5U8yLy8vJCcnP5avxDl69ChGjRqFixcvihbss/+RyWTIysr601/nwxgzHf6cKMa6obq6Gtu2bUNFRQVKS0sxf/58VFZW4qWXXjL11B7KkiVLYGNjg4aGhkd63qysLOTk5KCqqgqHDx/Ga6+9hpEjR3IDZUB8fDzf1mPMTPArUYx1w+XLlzFr1iycPXsWRISAgACkpqYKi7+loLq6WnhX2YABA4TF0o/Crl27sHLlSly+fBlOTk6YMGECNmzYAEdHx0d2jYc1ePDgTj+5fMuWLZ0uZv+r1dXVobGxEcDdj9Jo/45Fxpi0cBPFGDNL7ZvG+7m4uHT47CnGGHtY3EQxxhhjjBmB10QxxhhjjBmBmyjGGGOMMSNwE8UYY4wxZgRuohhjjDHGjMBNFGOMMcaYEbiJYowxxhgzAjdRjDHGGGNG4CaKMcYYY8wI/wfuNnLM7uctjwAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "bias = sample['tas'] - sample['2t']\n", + "bias.plot()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "92785d26-738a-4ccf-a0fa-56b7838a1cf0", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.7" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/packages/data/src/pyearthtools/data/transforms/transform.py b/packages/data/src/pyearthtools/data/transforms/transform.py index 4841121a..70f4af1c 100644 --- a/packages/data/src/pyearthtools/data/transforms/transform.py +++ b/packages/data/src/pyearthtools/data/transforms/transform.py @@ -242,7 +242,7 @@ def __call__(self, dataset: XR_TYPES | tuple[XR_TYPES] | list[XR_TYPES] | dict[s # Do not try to transform null datasets if dataset is None: return None - + # Do not try to transform empty datasets if not len(dataset): return None diff --git a/packages/nci_site_archive/src/site_archive_nci/BARRA.py b/packages/nci_site_archive/src/site_archive_nci/BARRA.py index 5df76dc5..0c1b1eec 100644 --- a/packages/nci_site_archive/src/site_archive_nci/BARRA.py +++ b/packages/nci_site_archive/src/site_archive_nci/BARRA.py @@ -136,7 +136,9 @@ def __init__( preprocess = None if datatype == "analysis": - base_transform += pyearthtools.data.transforms.coordinates.Drop(["forecast_reference_time", "forecast_period"]) + base_transform += pyearthtools.data.transforms.coordinates.Drop( + ["forecast_reference_time", "forecast_period"] + ) preprocess = pyearthtools.data.transforms.dimensions.Expand("time", as_dataarray=True) self.pressure = pressure diff --git a/packages/nci_site_archive/src/site_archive_nci/_CMIP5.py b/packages/nci_site_archive/src/site_archive_nci/_CMIP5.py index 9ef1bdb1..0cbd0a92 100644 --- a/packages/nci_site_archive/src/site_archive_nci/_CMIP5.py +++ b/packages/nci_site_archive/src/site_archive_nci/_CMIP5.py @@ -45,17 +45,21 @@ from pyearthtools.data.operations.index_routines import _mf_series # Path schema for CMIP5 -NamedPath = namedtuple("NamedPath", [ "institute", "model", "scenario", "interval", "realm", "cat", "expid", "expver", "variable"]) +NamedPath = namedtuple( + "NamedPath", ["institute", "model", "scenario", "interval", "realm", "cat", "expid", "expver", "variable"] +) + class CMIP_Path: - ''' + """ Interpret the path elements of a file containing CMIP data according to the path schema - ''' + """ + def __init__(self, *, elements=None, filepath=None): - + if not elements: dirpath, filename = os.path.split(filepath) - elements = dirpath.split('/')[-9:] + elements = dirpath.split("/")[-9:] self.namedpath = NamedPath(*elements) @@ -64,38 +68,58 @@ def __init__(self, *, elements=None, filepath=None): def __getitem__(self, keyword): return self.namedpath.__getattribute__(keyword) - -INSTITUTIONS = ['BCC', 'BNU', 'CCCma', 'CMCC', 'CNRM-CERFACS', 'FIO', 'ICHEC', 'INM', 'INPE', 'IPSL', 'LASG-CESS', 'LASG-IAP', - 'MIROC', 'MOHC', 'MPI-M', 'MRI', 'NASA-GISS', 'NASA-GMAO', 'NCAR', 'NIMR-KMA', 'NOAA-GFDL', 'NOAA-NCEP', 'NSF-DOE-NCAR'] -UNDER_DEV_MSG = ''' +INSTITUTIONS = [ + "BCC", + "BNU", + "CCCma", + "CMCC", + "CNRM-CERFACS", + "FIO", + "ICHEC", + "INM", + "INPE", + "IPSL", + "LASG-CESS", + "LASG-IAP", + "MIROC", + "MOHC", + "MPI-M", + "MRI", + "NASA-GISS", + "NASA-GMAO", + "NCAR", + "NIMR-KMA", + "NOAA-GFDL", + "NOAA-NCEP", + "NSF-DOE-NCAR", +] + +UNDER_DEV_MSG = """ Under development. This class currently only allows single values - e.g. one institution, one model, one experiment id etc. Expansion to allowing multiple values in under way, which will create a more complex DataSet object containing all relevant data in a single in-memory object. -''' +""" def rounder(time: Petdt, interval: int) -> str: hour = time.hour return "%02d00" % ((hour // interval) * interval,) - -SAMPLE_PLUGIN_DICT = { - 'models': 'bcc-csm1-1', - 'scenarios': "an", - 'institutions': 'BCC' -} + +SAMPLE_PLUGIN_DICT = {"models": "bcc-csm1-1", "scenarios": "an", "institutions": "BCC"} + def to_cftime(times): - ''' + """ Given an iterable of Petdts, return an iterable of cftime.DatetimeNoLeap - ''' + """ pdt = times[0] # converted = [cftime.num2date(pdt.datetime.timestamp(), 'seconds since 1970-01-01') for pdt in times] - converted = [pdt.to_cftime(calendar='noleap') for pdt in times] + converted = [pdt.to_cftime(calendar="noleap") for pdt in times] return converted @@ -145,28 +169,27 @@ def __init__( self.record_initialisation() - def quick_walk(self): - ''' + """ Walking a large filesystem to find matching filenames can take a long time. - This function uses the query dictionary to more effectively walk only + This function uses the query dictionary to more effectively walk only the parts of the filesystem actually relevant to the dataset. If performance is not a concern or if the filesystem is small, just use os.walk. - ''' + """ CMIP5_HOME = self.ROOT_DIRECTORIES["CMIP5"] basepath = Path(CMIP5_HOME) # Only walk the filesystem once, use the cache thereafter if self.walk_cache: - for (root, dirs, files) in self.walk_cache: + for root, dirs, files in self.walk_cache: yield (root, dirs, files) - return - + return + # Only walk the filesystem for configured institutions for institution in self.institutions: - institution_basepath = os.path.join(basepath, institution) + institution_basepath = os.path.join(basepath, institution) for model in self.models: model_basepath = os.path.join(institution_basepath, model) @@ -177,70 +200,65 @@ def quick_walk(self): all_entries = list(os.walk(scenario_basepath)) walk_cache = [] - for (root, dirs, files) in all_entries: + for root, dirs, files in all_entries: # We don't care about directories with no files if files: walk_cache.append((root, dirs, files)) - + self.walk_cache = walk_cache - def filesystem( - self, - query_dictionary={} - ): - ''' + def filesystem(self, query_dictionary={}): + """ Given the supplied query, return all filenames which contain the data necessary to extract the data for the query. For example, figure out what time index and variables the user wants, and go find all the files matching those time indexes and variables so they can get loaded into memory and then the relevant data extracted and re-aggregated as needed by other parts of the system. - ''' + """ paths = {} - paths_to_open = [] - + paths_to_open = [] # Walk the filesystem finding relevant file paths # A more efficient walk ordered by time or primary dimension may be more efficient - for (root, _dirs, files) in self.quick_walk(): + for root, _dirs, files in self.quick_walk(): paths = [os.path.join(root, file) for file in files] relevant = [p for p in paths if self.match_path(p, query_dictionary)] paths_to_open += relevant return paths_to_open - + def match_path(self, path, query): - ''' + """ Given a query (typically a date/time) and a full path to a filename, go and check if that filename is likely to contain data relevant to the query Returns: True if the path and query match False if the file should be ignored - ''' + """ path = CMIP_Path(filepath=path) match = True - if path['variable'] not in self.variables: + if path["variable"] not in self.variables: match = False - if path['model'] not in self.models: + if path["model"] not in self.models: match = False - if path['interval'] not in self.interval: + if path["interval"] not in self.interval: match = False - if path['scenario'] not in self.scenarios: + if path["scenario"] not in self.scenarios: match = False - if 'atmos' != path['realm']: + if "atmos" != path["realm"]: match = False # TODO: Filter by matching time range of filename or data - - return match + return match # Override the series method because of the unusual datetime class # Maybe the data loader should convert from unusual calendar datetimes to regular calendar datetimes @@ -260,7 +278,6 @@ def series( tolerance: tuple | pd.Timedelta | None = None, **kwargs, ) -> xr.Dataset: - """ Index into Provided Data function to create a continuous series of Data @@ -372,7 +389,9 @@ def series( and DataFunction.data_resolution > start.resolution ): timesteps = [ - t for time in timesteps for t in pyearthtools.data.TimeRange(time, time + 1, DataFunction.data_interval) + t + for time in timesteps + for t in pyearthtools.data.TimeRange(time, time + 1, DataFunction.data_interval) ] time = list(set(map(lambda x: x.datetime64("ns"), timesteps)) & set(data[time_dim].values)) @@ -434,4 +453,4 @@ def series( IndexWarning, ) - return transforms(data) \ No newline at end of file + return transforms(data) diff --git a/packages/nci_site_archive/src/site_archive_nci/__init__.py b/packages/nci_site_archive/src/site_archive_nci/__init__.py index 86867d83..692e4c2c 100644 --- a/packages/nci_site_archive/src/site_archive_nci/__init__.py +++ b/packages/nci_site_archive/src/site_archive_nci/__init__.py @@ -35,7 +35,7 @@ from pyearthtools.data.archive import register_archive # Please note, these directories may be found on the public NCI documentation -# site, in the geonetwork data catalogue. +# site, in the geonetwork data catalogue. # See e.g. https://geonetwork.nci.org.au/geonetwork/srv/eng/catalog.search#/metadata/f7332_1490_5979_5715 ROOT_DIRECTORIES = { @@ -50,7 +50,7 @@ "BARRA": "/g/data/cj37/BARRA/BARRA_{region}/{version}/{datatype}", "BARPA": "/g/data/py18/BARPA/", "BARRA_V2": "/g/data/ob53/BARRA2/", - "CMIP5" : "/g/data/al33/replicas/CMIP5/combined/", + "CMIP5": "/g/data/al33/replicas/CMIP5/combined/", } diff --git a/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/__init__.py b/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/__init__.py index 6ba4a0f4..c120c85f 100644 --- a/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/__init__.py +++ b/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/__init__.py @@ -37,6 +37,8 @@ from pyearthtools.pipeline.operations.xarray.join import Merge, Concatenate from pyearthtools.pipeline.operations.xarray.sort import Sort from pyearthtools.pipeline.operations.xarray.chunk import Chunk +from pyearthtools.pipeline.operations.xarray._recode_calendar import RecodeCalendar +from pyearthtools.pipeline.operations.xarray._align_dates import AlignDates from pyearthtools.pipeline.operations.xarray import ( conversion, @@ -65,4 +67,6 @@ "metadata", "normalisation", "remapping", + "RecodeCalendar", + "AlignDates", ] diff --git a/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/_align_dates.py b/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/_align_dates.py new file mode 100644 index 00000000..c262e4a5 --- /dev/null +++ b/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/_align_dates.py @@ -0,0 +1,77 @@ +# Copyright Commonwealth of Australia, Bureau of Meteorology 2025. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + + +from typing import TypeVar, Optional + +import pandas as pd +import xarray as xr + +from pyearthtools.pipeline.operation import Operation + +T = TypeVar("T", xr.Dataset, xr.DataArray) + + +class AlignDates(Operation): + """ + Climate datasets often use the cftime module to index into data using non-standard calendars. + This operation will recode the time coordinate of a dataset or data array to a standard timestamp + For now, support only exists for recoding from Noleap to Timestamp + """ + + _override_interface = "Serial" + + def __init__(self, to="start_of_month"): + """ + Record initialisation and store flags for processing + + Args: + to: either "start_of_month" or a zero-padded string day-of-month + """ + + if to == "start_of_month": + to = "01" + + if len(to) != 2: + raise ValueError(f"Value of 'to' {to} is not recognised as a day of the month") + self.to = to + + super().__init__( + split_tuples=True, + operation="apply", + recognised_types=(xr.Dataset, xr.DataArray), + ) + self.to = to + self.record_initialisation() + + def apply_func(self, data: xr.Dataset) -> xr.Dataset: + """Sort an `xarray` object data variables into the given order + + Args: + data (T): + `xarray` object to sort. + + Returns: + (T): + Sorted dataset + """ + to = self.to + if self.to == "start_of_month": + to = "01" + + aligned = data.time.dt.strftime(f"%Y-%m-{to}") + aligned = [pd.Timestamp(a) for a in aligned.values] + data = data.assign_coords({"time": aligned}) + + return data diff --git a/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/_recode_calendar.py b/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/_recode_calendar.py new file mode 100644 index 00000000..c37f83bb --- /dev/null +++ b/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/_recode_calendar.py @@ -0,0 +1,61 @@ +# Copyright Commonwealth of Australia, Bureau of Meteorology 2025. +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. + + +from typing import TypeVar, Optional + +import xarray as xr + +from pyearthtools.pipeline.operation import Operation + +T = TypeVar("T", xr.Dataset, xr.DataArray) + + +class RecodeCalendar(Operation): + """ + Climate datasets often use the cftime module to index into data using non-standard calendars. + This operation will recode the time coordinate of a dataset or data array to a standard timestamp + For now, support only exists for recoding from Noleap to Timestamp + """ + + _override_interface = "Serial" + + def __init__(self): + """ + Record initialisation and store flags for processing + """ + + super().__init__( + split_tuples=True, + operation="apply", + recognised_types=(xr.Dataset, xr.DataArray), + ) + self.record_initialisation() + + def apply_func(self, data: xr.Dataset) -> xr.Dataset: + """Sort an `xarray` object data variables into the given order + + Args: + data (T): + `xarray` object to sort. + + Returns: + (T): + Sorted dataset + """ + + recoded = data.indexes["time"].to_datetimeindex() + data["time"] = recoded + + return data diff --git a/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/join.py b/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/join.py index eec4eb0f..6a527438 100644 --- a/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/join.py +++ b/packages/pipeline/src/pyearthtools/pipeline/operations/xarray/join.py @@ -45,6 +45,68 @@ def unjoin(self, sample: Any) -> tuple: return super().unjoin(sample) +class GeospatialTimeSeriesMerge(Joiner): + """ + The default "merge" and "interplike" xarray commands are very general + This can result in 'merged' or 'interpolated' arrays which don't do + what's hoped for but don't raise any errors either. + + This joiner is more strict about the merging and interpolating, and also + raises more informative error messages when it runs into trouble. + """ + + _override_interface = "Serial" + + def __init__( + self, + reference_dataset=None, + reference_index=None, + interpolation_method="nearest", + time_dimension="time", + merge_kwargs: Optional[dict[str, Any]] = None, + ): + super().__init__() + self.record_initialisation() + self.reference_dataset = reference_dataset + self.reference_index = reference_index + self.interpolation_method = interpolation_method + self.time_dimension = time_dimension + self._merge_kwargs = merge_kwargs + + def _join_two_datasets(self, sample_a: xr.Dataset, sample_b: xr.Dataset) -> xr.Dataset: + """ + Used to reduce a sequence of joinable items. Only called by the public interface join method. + """ + + # Check each sample has the proper time dimension + if self.time_dimension not in sample_a.coords: + raise ValueError(f"Time dimension missing from {str(sample_a)}") + + if self.time_dimension not in sample_b.coords: + raise ValueError(f"Time dimension missing from {str(sample_b)}") + + interped_a = sample_a.interp_like(self.reference_dataset, method="nearest") + interped_b = sample_b.interp_like(self.reference_dataset, method="nearest") + merged = xr.merge([interped_a, interped_b]) + return merged + + def join(self, sample: tuple[Union[xr.Dataset, xr.DataArray], ...]) -> xr.Dataset: + """Join sample""" + + # Obtain the reference dataset + if self.reference_dataset is None: + if self.reference_index is not None: + self.reference_dataset = sample[self.reference_index] + else: + raise ValueError("No reference dataset or reference index set") + + merged = reduce(lambda a, b: self._join_two_datasets(a, b), sample) + return merged + + def unjoin(self, sample: Any) -> tuple: + raise NotImplementedError("Not Implemented") + + class InterpLike(Joiner): """ Merge a tuple of xarray object's. @@ -54,17 +116,32 @@ class InterpLike(Joiner): _override_interface = "Serial" - def __init__(self, reference_dataset=None, merge_kwargs: Optional[dict[str, Any]] = None): + def __init__( + self, + reference_dataset=None, + reference_index=None, + method="nearest", + merge_kwargs: Optional[dict[str, Any]] = None, + ): super().__init__() self.record_initialisation() self.reference_dataset = reference_dataset + self.reference_index = reference_index + self.interp_method = method self._merge_kwargs = merge_kwargs def join(self, sample: tuple[Union[xr.Dataset, xr.DataArray], ...]) -> xr.Dataset: """Join sample""" # merged = reduce(lambda a, b: a.interp_like(b), sample) - reference = self.reference_dataset - interped = [i.interp_like(reference) for i in sample] + + if self.reference_dataset is not None: + reference = self.reference_dataset + elif self.reference_index is not None: + reference = sample[self.reference_index] + else: + raise ValueError("No reference dataset or reference index set") + + interped = [i.interp_like(reference, method=self.interp_method) for i in sample] merged = xr.merge(interped) return merged