{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Calculate Time Averages from Time Series Data\n", "=============================================\n", "\n", "Author: [Tom Vo](https://github.com/tomvothecoder/)\n", "\n", "Date: 05/27/22\n", "\n", "Last Edited: 08/17/22 (v0.3.1)\n", "\n", "Related APIs:\n", "\n", "* [xarray.Dataset.temporal.average()](../generated/xarray.Dataset.temporal.average.rst)\n", "* [xarray.Dataset.temporal.group_average()](../generated/xarray.Dataset.temporal.group_average.rst)\n", "\n", "The data used in this example can be found through the [Earth System Grid Federation (ESGF) search portal](https://aims2.llnl.gov/metagrid/search)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Overview\n", "\n", "Suppose we have netCDF4 files for air temperature data (`tas`) with monthly, daily, and 3hr frequencies.\n", "\n", "We want to calculate averages using these files with the time dimension removed (a single time snapshot), and averages by time group (yearly, seasonal, and daily)." ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "ExecuteTime": { "end_time": "2018-11-28T20:51:35.958210Z", "start_time": "2018-11-28T20:51:35.936966Z" } }, "outputs": [], "source": [ "%matplotlib inline\n", "\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "import xcdat\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Calculate averages with the time dimension removed (single snapshot)\n", "\n", "Related API: [xarray.Dataset.temporal.average()](../generated/xarray.Dataset.temporal.average.rst)\n", "\n", "Helpful knowledge:\n", "\n", "* The frequency for the time interval is inferred before calculating weights.\n", " * The frequency is inferred by calculating the minimum delta between time coordinates and using the conditional logic below. This frequency is used to calculate weights.\n", "\n", " ```python\n", " if min_delta < pd.Timedelta(days=1):\n", " return \"hour\"\n", " elif min_delta >= pd.Timedelta(days=1) and min_delta < pd.Timedelta(days=28):\n", " return \"day\"\n", " elif min_delta >= pd.Timedelta(days=28) and min_delta < pd.Timedelta(days=365):\n", " return \"month\"\n", " else:\n", " return \"year\"\n", " ```\n", "* Masked (missing) data is automatically handled.\n", " * The weight of masked (missing) data are excluded when averages are calculated. This is the same as giving them a weight of 0." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Open the ``Dataset``\n", "\n", "In this example, we will be calculating the time weighted averages with the time dimension removed (single snapshot) for monthly `tas` data.\n", "\n", "We are using xarray's OPeNDAP support to read a netCDF4 dataset file directly from its source. The data is not loaded over the network until we perform operations on it (e.g., temperature unit adjustment).\n", "\n", "*More information on the xarray's OPeNDAP support can be found [here](https://docs.xarray.dev/en/stable/user-guide/io.html#opendap).*" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
       "Dimensions:    (time: 1980, bnds: 2, lat: 145, lon: 192)\n",
       "Coordinates:\n",
       "  * time       (time) datetime64[ns] 1850-01-16T12:00:00 ... 2014-12-16T12:00:00\n",
       "  * lat        (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n",
       "  * lon        (lon) float64 0.0 1.875 3.75 5.625 ... 352.5 354.4 356.2 358.1\n",
       "    height     float64 2.0\n",
       "Dimensions without coordinates: bnds\n",
       "Data variables:\n",
       "    time_bnds  (time, bnds) datetime64[ns] 1850-01-01 1850-02-01 ... 2015-01-01\n",
       "    lat_bnds   (lat, bnds) float64 -90.0 -89.38 -89.38 ... 89.38 89.38 90.0\n",
       "    lon_bnds   (lon, bnds) float64 -0.9375 0.9375 0.9375 ... 357.2 357.2 359.1\n",
       "    tas        (time, lat, lon) float32 -27.19 -27.19 -27.19 ... -25.29 -25.29\n",
       "Attributes: (12/49)\n",
       "    Conventions:                     CF-1.7 CMIP-6.2\n",
       "    activity_id:                     CMIP\n",
       "    branch_method:                   standard\n",
       "    branch_time_in_child:            0.0\n",
       "    branch_time_in_parent:           87658.0\n",
       "    creation_date:                   2020-06-05T04:06:11Z\n",
       "    ...                              ...\n",
       "    version:                         v20200605\n",
       "    license:                         CMIP6 model data produced by CSIRO is li...\n",
       "    cmor_version:                    3.4.0\n",
       "    _NCProperties:                   version=2,netcdf=4.6.2,hdf5=1.10.5\n",
       "    tracking_id:                     hdl:21.14100/af78ae5e-f3a6-4e99-8cfe-5f2...\n",
       "    DODS_EXTRA.Unlimited_Dimension:  time
" ], "text/plain": [ "\n", "Dimensions: (time: 1980, bnds: 2, lat: 145, lon: 192)\n", "Coordinates:\n", " * time (time) datetime64[ns] 1850-01-16T12:00:00 ... 2014-12-16T12:00:00\n", " * lat (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n", " * lon (lon) float64 0.0 1.875 3.75 5.625 ... 352.5 354.4 356.2 358.1\n", " height float64 2.0\n", "Dimensions without coordinates: bnds\n", "Data variables:\n", " time_bnds (time, bnds) datetime64[ns] ...\n", " lat_bnds (lat, bnds) float64 ...\n", " lon_bnds (lon, bnds) float64 ...\n", " tas (time, lat, lon) float32 -27.19 -27.19 -27.19 ... -25.29 -25.29\n", "Attributes: (12/49)\n", " Conventions: CF-1.7 CMIP-6.2\n", " activity_id: CMIP\n", " branch_method: standard\n", " branch_time_in_child: 0.0\n", " branch_time_in_parent: 87658.0\n", " creation_date: 2020-06-05T04:06:11Z\n", " ... ...\n", " version: v20200605\n", " license: CMIP6 model data produced by CSIRO is li...\n", " cmor_version: 3.4.0\n", " _NCProperties: version=2,netcdf=4.6.2,hdf5=1.10.5\n", " tracking_id: hdl:21.14100/af78ae5e-f3a6-4e99-8cfe-5f2...\n", " DODS_EXTRA.Unlimited_Dimension: time" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "filepath = \"http://esgf.nci.org.au/thredds/dodsC/master/CMIP6/CMIP/CSIRO/ACCESS-ESM1-5/historical/r10i1p1f1/Amon/tas/gn/v20200605/tas_Amon_ACCESS-ESM1-5_historical_r10i1p1f1_gn_185001-201412.nc\"\n", "\n", "ds = xcdat.open_dataset(filepath)\n", "\n", "# Unit adjust (-273.15, K to C)\n", "ds[\"tas\"] = ds.tas - 273.15\n", "\n", "ds\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "ds_avg = ds.temporal.average(\"tas\", weighted=True)" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'tas' (lat: 145, lon: 192)>\n",
       "array([[-48.01481628, -48.01481628, -48.01481628, ..., -48.01481628,\n",
       "        -48.01481628, -48.01481628],\n",
       "       [-44.94085363, -44.97948214, -45.01815398, ..., -44.82408252,\n",
       "        -44.86273067, -44.9009281 ],\n",
       "       [-44.11875274, -44.23060624, -44.33960158, ..., -43.76766492,\n",
       "        -43.88593717, -44.00303006],\n",
       "       ...,\n",
       "       [-18.21076615, -18.17513373, -18.13957458, ..., -18.32720478,\n",
       "        -18.28428828, -18.2486193 ],\n",
       "       [-18.50778243, -18.49301854, -18.47902819, ..., -18.55410851,\n",
       "        -18.5406963 , -18.52413098],\n",
       "       [-19.07366375, -19.07366375, -19.07366375, ..., -19.07366375,\n",
       "        -19.07366375, -19.07366375]])\n",
       "Coordinates:\n",
       "  * lat      (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n",
       "  * lon      (lon) float64 0.0 1.875 3.75 5.625 7.5 ... 352.5 354.4 356.2 358.1\n",
       "    height   float64 2.0\n",
       "Attributes:\n",
       "    operation:  temporal_avg\n",
       "    mode:       average\n",
       "    freq:       month\n",
       "    weighted:   True
" ], "text/plain": [ "\n", "array([[-48.01481628, -48.01481628, -48.01481628, ..., -48.01481628,\n", " -48.01481628, -48.01481628],\n", " [-44.94085363, -44.97948214, -45.01815398, ..., -44.82408252,\n", " -44.86273067, -44.9009281 ],\n", " [-44.11875274, -44.23060624, -44.33960158, ..., -43.76766492,\n", " -43.88593717, -44.00303006],\n", " ...,\n", " [-18.21076615, -18.17513373, -18.13957458, ..., -18.32720478,\n", " -18.28428828, -18.2486193 ],\n", " [-18.50778243, -18.49301854, -18.47902819, ..., -18.55410851,\n", " -18.5406963 , -18.52413098],\n", " [-19.07366375, -19.07366375, -19.07366375, ..., -19.07366375,\n", " -19.07366375, -19.07366375]])\n", "Coordinates:\n", " * lat (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n", " * lon (lon) float64 0.0 1.875 3.75 5.625 7.5 ... 352.5 354.4 356.2 358.1\n", " height float64 2.0\n", "Attributes:\n", " operation: temporal_avg\n", " mode: average\n", " freq: month\n", " weighted: True" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds_avg.tas" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEXCAYAAABcRGizAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAACVEklEQVR4nO39eZQsy1Uein87MrOquvuce+4odDVZVzKDEWBJCDE9Y2EGCzEPFmCDmdYTg7EBG7CEsI2RwQgxmYfBXGz89AwIsEEgMCAk/0DMgyQkISEJhLhoutz5TN1dVZkZ+/dHDLkjcqis6uruc/rEt1avrpwiIjMjY8f+9hDEzEhISEhISBiCOu0GJCQkJCRc+0jCIiEhISFhJZKwSEhISEhYiSQsEhISEhJWIgmLhISEhISVSMIiISEhIWElkrBI6AUR3UNEn7jhtW8homcddz3XKoiIiWifiL5jw+s/kYiuEpE+a88m4fpEEhYJxwJmfgoz/+ZRyyGiZxHRe7bQpK6yv4eI/oKIrhDR24jon644/x8T0V9bIfALRHTriir+LjO/cJO2MfOrmfkcgHdtcn1CwraRhEXCjYx9AJ8O4AKALwHwn4joY7pOJKKnAPhRAF8M4P0AHAD44RNqZ0LCqSMJi4RVeCoRvYmILhHRzxDRzB0gok8jojcQ0UUi+j0i+jBxzFNLRLRDRC8lokeI6K1E9M0d2kKrHiLaA/CrAB5jKZmrRPSYbd0YM/87Zn4bM2tm/kMAvw3go3tO/ycAfomZf4uZrwL4NwA+h4jOj6mLiJ5oqakvI6J322fxVUT0Efa+LxLRD23nzhISto8kLBJW4bkAng3gLgAfBuBLAYCIng7gxwF8JYDbYGbdryCiaUcZ/w7AEwE8CcAnAfiiMfUw8z6ATwHwPmY+Z//eF19IRM+3g23n35ibJKIdAB8B4C09pzwFwBvdBjP/JYAlgA8YU77ARwJ4fwCfD+AHALwQwCfa8p9LRH9/zfISEk4ESVgkrMIPMvP7mPlhAL8E4Kl2//8N4EeZ+Q+ZuWbmlwJYAPiojjKeC+A7mfkRZn4PgB9co56VYObvYuab+/5GFvNfYITBK3uOnwNwKdp3CcAozULgRcw8Z+Zfh6HBXsbM9zPze2E0m6etWV5CwokgCYuEVfgb8fsAZtAEgL8F4F9FM/jHA+iiiR4D4N1i+90d5/TVc+wgopcA+BAAz+X+zJpXAdwU7bsJwJU1q7tP/D7s2D6x+05IWAdJWCRsincD+I5oFr/LzC/rOPdeAI8T249fo56VaZGJ6FuETaP1t+Lafw9DdX0yM18eOPUtAP6uuO5JAKYA/nzcbSQkXN9IwiJhU/wYgK8ioo8kgz0i+tQeg+/PAngBEd1CRI8F8LVr1HMfgNuI6ELfCcz8ncKm0frru46IXgDgHwP4JGZ+aEU7fhLApxPR37OG928H8PPMvK5mkZBwXSIJi4SNwMyvhbFb/BCARwC8A9b43YFvB/AeAH8F4NUA/heMfWNMPW8D8DIA77R019a8oQB8J4AnAPgLoYl8iztot/+ebcdbAHwVjNC4H8ZW8TVbbEtCwjUNSosfJZw0iOirAXwBM59Zzx8imsMIxB9k5n+zwfWfAODnYKiu5zDzb2y5iQkJayEJi4RjBxHdCeM2+/swbqP/G8APMfMPnGa7EhISxiM/7QYk3BCYwMRh3AXgIoCfRop+Tki4rpA0i4SEhISElUgG7oSEhISElTgzNNS5m2/l2+583OoTExISbni8621/+iAz33GUMh5POzyHHnXug1i+kpmffZT6ThtnRljcdufj8Pz//orTbkZCQsJ1gK/56Lv++qhlLKDxXLpz1Lk/zH99+1HrO22cGWEBAIqoc3/WvTthJOrrxKylk/3tukffN3wtggBkY9t7BrrmmREW5sV1H1Pq+umApwWt+3vzusL2OITLSQqCOgmdrWPsoNr1njcVICcxSbyRJqJnRlgkJCQknCTW0izOAM6OsKBGgxj7Ao9D5R2aaZwUnbPJLDwTDV93Zh1rJWNmW+s+C/eutqVhnJT2UA9obDcSasHDZCM0ffkNx+98E7r5WNgFSprFdQkCUKi2J3DXyxzbcbYtTMY87E0HQzlgDwnLMYO0IlrZDjnY9j3PMdTWJkLjKAJjW0LitITAtV6vEwS15l6hsM49yDJ82fYdyn6uiFBz+L3H/XLbWoACMLmBKO4zIywSEhISThaUaKjrEc7A7WYTUivoNXwPvOgxqvI20J5l9czSV830R3RazdwbhdmmkiQttbLosC22rV0axxjKalV9R9UujoKTmNlvs47jbm9cvtweqrvWPKptsWYxyVXvMdNnHRVt90fahzx2VBBurKjmsyMsiJBnyggMoqAjDQmOMVTTtgVH8JFk5t9q2qf/mGZGsaoeAED/IBsLG3me/CDcYO8+wrhdmrk12+qjrPpoqozGCYy4ndczjmtQ76RxYoFt9y8r3XkdAD9Iy2uXlUatORjA4zJcObK8Sa5wuKw773loX60Zy0r7snYmWdD2QJBEfbBQ1DmRPCqSZpGQkJCQMAhKBu7rE4qA3SLzmkWRNZ5RXv1U/bOKVTPUfg+M4dlKXG6pNcqBafOQAXbdyWetGZo5mKVLw3KgPUSFx3W5MlQWzugV2sb1MbRVzTyoZcQfYV+Z8rmP0TIajehsaCRDkEZm+TuehbtjO5PMH3e4aZpjv6xxuKwBINAIHC3ktg+Xta/jwm6o6y4rHZTvfgOrHDLatFaX9rFTZCgkm2BpKbfPsQ4AUGTbIY8IycB9XSIjwoVpjiIj7BYZcvsSC0VwfYMAHFZhR1vW5sM4KOtg/zTL/HVuQJJ9WgnbQl9fr5nBHB5UNQGouy8AUHTYLGwTBwfDLg8RACgrRmkLKDLl76XIKPjotAqFSstdUZ4rqKaajQAJzh+ww+ionX32jSHhMUZwdN1DWN5mrsJ9VM5RIWmadcvuo24yRbg6r3z5TkDcvjfxgylg7l9rxiOHZauc/dIIh8feNANgvpPS1nfpsMSFnQKX7HXnZs1wEntD7UyyTgqr6/4D+kpSyLYT1Kq5znlAFllobC6UoaWdsCiOSVgkGiohISEhYSUSDXUdIleEO/YKTDPCJFNeKzA0lPmtGTgoKz+TXFTa/94tMtQanr7KhUYiywD6/JUauHmPZgIzoO0eZlPWJCe/7Wa/dZS8Us6KC+Vmv91G3Vp302CFAlShzA97/242V9Y6cAQoiLxCYM4JtQVpUM7QaCWKujUNeR9uHqc19xrGg/uJKCp3rUPfBxqXOVbT6JsdDmkcQ04PR9U61tEy4nbE3kiO7skU4dKB0QAm51Wg/dXMePjq0m9PcoXHXtgBAJybZDgsa/9sp3mGcmm0lVt2J9DMuLBT+HIkVXS4rD2ltay0/72wWoU0fjsmYGeSYZIrrwVNcoWJ/RDj91kowjQ391dqHdDP7lijWRBmtszdIsM2YGwWN460ODPCIlOEC1NDP2UEocaG5z32fI7LS0c9EXbtQKoAlJr9QD8RqnDcH4ZYC82AW1BKw3w07Ad50x4nF8xpzTFZds0cCBLpLhp3UEV9Der3fiprtnU0+/yHpcgeExywPdGd32cr6KKrpL0jtnVsKjx8HSOEiCt/rBdM33Meap/Ecbtdx3RO17EuLyRHEx2WbRp0WWlvZ3jsTTPcNDO/MwJmucK7Ls39thtsDX0FT0vJ9iwrjavzCgfW1nF1UeHQCpmFEByujbtWqO1Mcpyb5jhv23puliObNXZBIwyae1pUpvxpnqHICFMrWGa5wjRXgYBwQmcn36Y31NaKuuZxqsKCiD4QwM+IXU8C8G8B3Azg/wbwgN3/Lcz8KyfbuoSEhIR+GJvFabfi5HCqwoKZ3w7gqQBARBmA9wJ4OYAvA/D9zPw9Y8tSZGYMLjjHTQopnvXlCrfa2dZNE9XrYaQonJVJjcFt+/uA1AwYGs4ICijVzJoJDCaA7Iy9Bgf007LWfrvIGi8uWINe1kOAKS2Ng83+WsN6QzWNdZ5YZRQUVQMoo/w9wQTZzVQtBdVQe/3aQ0GqFQPSZ1RX1sgKtD2qBoPwxtoq10kxMULbCZpwxJgPRUYbdoZX+b4WVcfiOpZFWdYm5mAiXvoYY31XQNyjb57hwtS8KaXIz9gLayS+8/wUAHBpXnmqx8H1qYOyRjY1xw5VGEexrDWWVeONdbisUdt7U4q8t9W5mdk37YiZmOUq8HiqGcgEo1SokGo6N8mxUzhNo9E6Jlsa4Qm0VW8oOwa+FsB7mfnTiOhWmMn0EwHcA+C5zPzI1ipcE9cSDfUJAP6Smf+aNuABCeajUWCAa5C2Hxnr6DxAkek0uVKA/c1EoWBhDW99IAVW5D2bNHOQnl4KkpoJnhVS1hvKfjREhMNK47Bs2hQEDAo7iUQXHSKbmuXNdlkzSnvv7iP2dNagW254TNfsP1hTfzB6ixN9i8TVgr4aGOhVZE+qxaDb5ZLrEAuSofvw128wQDghN/jcRubkiiFpNeneDUR9ojCUYFfZXdy7bGtf0GOtwudUM+P8JPdCQGtG6Y9plNS06cIsbzzqFAW2gEvzyntU7RYZdorMU16TXGGnaKimw2UV2DcaG4W5TtJorr5ph7Bw25oZ0zzzz+TW3QI7ucJO3giIwtkgtzTAH4PN4usAvBXATXb7+QD+DzN/FxE9327/621WuA6upWj1LwDwMrH9tUT0JiL6cSK6pesCInoeEb2WiF774IMPdJ2SkJCQcCxwNNSYv5VlET0OwKcC+K9i92cCeKn9/VIAn7XdO1gP14RmQUQTAJ8B4AV2148AeBHMdPVFAL4XwJfH1zHz3QDuBoAPf/rTWXEN0lU47WYN4g5VHuZls9UsWu+ThBwlApHy5ypSYGo8mgy91EyRm7IIDPYz8bLSWNYsvDvGsSOd55F7Bsbbys0U5YzRaylazsZCqqnxagqfQDz7ioMPPdVkjeHSUN8XkxGjlXbCn2br6vjKhvJbuWvljHpMhuFVhm/dMVv32k1HY8Ys4CPzmLnZsTQUu3rawZHjtJwhCq3UoTcSAMwrHcTjOK1Kg6EVBbRjcy8quNcLs9xrow8elCbFhqP0dOgpBbSN80DbMJ+pxsPJGLibB76ryMdHTfMMF6Y5brWeWecmCtO8oYlkHIba4rJ1W9QsfgDANwM4L/a9HzPfCwDMfC8RPWpblW2Ca0JYAPgUAK9n5vsAwP0HACL6MQC/vLoIBtVlIxykgOBuSgqkQiERCIhGOIABEIPIXM8qB7lRgghKeB0RkRdWisxHW4sBZpapwLvKfWs66sDGrdYe45D5cftMmTZK295amAfLDUbNuc315hqt2h9OPCB1DVCl+PANBWCDo2JbRxCQZ2gN07bIwytyx+3DKNNDT14wIBSCm6yLEEe+u/uXx/sQpM+OBr5FVQcUi7uHnLuFlazP0Y7OjRToFljuPRZZSF/FA96VZeXpnFnm3K7bHmJuIiCb5GmgHeDhwxIukFtHbrVdwgNoggel66wPvIvdhLm5Zycozk+tjSIjTLLGJTdTBNJGsPRNHteF0yxG4nYieq3YvttOdkFEnwbgfmZ+HRE9ayuNOwZcK8LiCyEoKCK600lUAJ8N4M2n0qqEhISEAayhWTzIzM/oOfaxAD6DiJ4DYAbgJiL6CQD3ubGQiO4EcP/RW7w5Tl1YENEugE8C8JVi93cT0VNh5vT3RMe6wWwpKN38Oei+mYR1V3JQeUNL6QpQ3Y+HWOgB5HzOXTNYBOWZ/242O0M489RsvKMAM7OWM7SKOZhhB95X0rhtPb8C7yRH43QYf2XqEBOfEc4UAZNyxHlSxeU47yrfNkWBl1aXMdbNaGvNw1l/Vce+NRHP3lvHA81CXjescYxJtWLOg6dzzLGwXbFDg8MsytXk6is7VKnaU1TWm8yEL0ReSh33MEBTSe+7cyJvU6nbfcDB7ZPP0d2HaX+BBw9MsN80Vz7Ow2sW3KajJpnC7iTzMR8XpoWPhYpTrdRgnJ+aMi/McuwW5D2dchsn5WOltKWoAUD3p9tZB0SELD+62ZeZXwBLwVvN4huZ+YuI6CUAvgTAd9n/v3jkyo6AUxcWzHwA4LZo3xdvUBLghIWO7BQd9BMAcJYD2cT+LlCqCRb2K5pOCJXtlI2XlW6XwS4HVDfHDDQDBWWETNBL5lrY67h38NJoD7L+1ti44nLHde4nc3PMBfAZ28KQSyp7W4dMwKgZUWqrkJ6ohatu/CwyRcJ7qotuagu3PvQJgpgWkgNoLBBk7q/YI6nUOrT/CAEcJ6d0qLW5blE1NojpBgLJPBdbX82t5xAL71jQxOX5Y+K6vON4VzllrQPhoVlEXlPzPFy73LHdIkNZMx61Z76v+/ebCPHlJLSZLCrt6aJzsxwXdgvcumOuuzDLvcsrENrZTCCuGcJ2bKS3jNjOJfWkKy8ktkVDgbrtalvEdwH4WSL6CgDvAvCPjrOyVTh1YZGQkJBwPYIA0JaSEjow828C+E37+yGYkIJrAmdHWDADddUYt3tnzBl4sgsAqKfnmpgHBrhuPJmWdTPTd5lpXfqLIhOxFIhnfTIlBnyAHmADz6hxoJHGUkcI+QA+OSMdyEbFAIi4MZRHdFVsOPflspkVu7JrDaEROEOom0ESQi8oHRjU5feiNHmDa5ehum8m1qmR+baFBtC47GBmG60CJYMQZQZiReQNp12p6zOlgn0+PhKmjC6mrNbAhJtypXYiZ+B9kM/fPd8sSuUSU16mXslLumO8Mg291ORkVlYJF1chHREKJZ53jzfCLFeoJzkWtZnN7xaZp8FcvIVMk+5yWJ2f5rgwzb2hXKbskO/MlKN86o6dQmEigvKUc3hxlBPrhobaomZBN1AI99kRFkDjNis6g1Q5mRQ43wHnJhpVCgR/fgcPCzR5nYDxieLGdCP5sZnvUXzcA4OLHwQIYGHvqLSwS4CDZIVxgJwSke4mYt2thha3UVZsPliZj8pEXDc8cpz+u2/Nj4DmoX63UBWln44pJNmOdrI59u0IkzOKZJOKEPnFBYkkFag3P5h7xr5tIEzsV1VwFjz73jKErQhAy0uty4tNbrvBOzyPrHCg1nlNuc1veazICFNx3kGpg+Nza2OY5QrQ8DRRLAw1M6bW80pRFaxUN8kV9uyDMtHWmf/t8jq57YmYtGSqWUNikoU2ipwAqkwOK0M7VU1jRtkw1wTRcdNQ1xTOlLBISEhIOCkQAdmWMtheDzg7woIZ0HU7xkJC5eB8hpLMbSsg4BMKQnNtR+Cdh7gmngtLDYSt0dorImRKcZoJZRTMEjmalfelPeHI40ijmaUqAko7669qQimeRRhQZW7Dz6ipmRlmUeyFnHVnys5sI4cS354MiJOahh5AzT6poMVp1vuM285rR2orYwzh7jy5+JOklqRHGRA+m5iW1tx0gQwUqJDSacGW1Nu2xtEhpsAajcCV2VwTerRJjcl5sbljTcBkux0yeNKn+ZBslngXs8jjx8WHHJQ1bt/Ng2fqrptkhL1J5ilcGaCXEZBlCuescfqWWYHzU0s7ZSZ7rEvTkSkEC5nlIo5H3hHDpOxRTpso5+0A3WNAoqGuSwjX2RguuE5l4HzqOzSLj15VC8NvinPdbyLVGyjGREF/lELFCRvyg4JJCe7d+RhQwcBMXkAESQvjOqPBhdgIGtNW8h9RRkDOqtOeAfQHvxG3DygxOJpBNMw/FdA/YrI1xNhJ/llRSBPpuikzU825sY1ErkFS63BgD2i2DipJrleSUTu7lRzchhC/Hzdgh0N+WxhydH5zmFo0VxOEKd8b9bhDC7uZ8GiT52VCOMb0nfQga3lUMft3ccusCJ53oVRAC+3kJldaXOY0zzDNlY+2ftTexAsHc20joKZZaIeQGRmW1AxfZc2YcEM9UV02Ngpge3YKCaKtG7ivZZwhYZGQkJBwciAcu+vsNYWzIyyYm1kF6yZVh8rBhVn1q967zc5M7CHWoGphNnRls8vaR6Lh03sEaUBiuJQhTgtR1KINPC1FxovKaQ3mOmm4ldc1Okq8thGDAw3ELKjUaEuutRnMbNJd3hfYFyPu/60ZuTCGG1qqO91IbWmfIS8gnztIUERlbX772WUepkgJvaPQ+u3bkiGgRfo0qVyZ5C1hzEuTKiKmp4AwJZRczEpzGOfpjikYrS/UCmx77TOUa70378ctRGXPZQSaIjN52pC5oeRcUGUh2nJuonz5sZOGpEO7UpA7r6YCjRbiVnxs0so0nlKFIswygou3y8T7zAg4N8mDGAmnWMxyo514Gqo8AC0iLyb7rU3FdzllDXX1UtPoruDcbYPM936j4MwICwJAte1USoEL4x6rdy547yfAeke54JzyoLkGJkiv+ZZFbijUiBMLBuWRCjoyiXNlZ3K2D0dPDZnGKBjYQ2ohsJlEZUk+3a2nEa7MNx5dY7zzrnJUVUaEWrHg3xt33VqzH/iB9iBk2uvureHXi8x4XLnzA9sK2unGZZpzANibNYNiE5DYf5+ZMhRV8zFQKz+g9JKLi5LCmJk99aK5eS9dY4qnwUBBwFzbE4/E6otNUBrbOpyQq3VDUWY2Pb6k4YK6s3BSUygEy/9K7E0UdrgJQs2jjiuFhROwROaZuVxNt1RF4H4tl0Cd5kawANYFlivQ4qopvFoKl1c7yZIusFHet8GAu215QXkQVKKhEhISEhKGQApQkyQsrkNI/VyDC6tNiCgtp1V4X+xqGU6jSMEn2unKSOsXSgqz04JrPw33moY7LM4nUi2aKqgimgk3sRSGkpLaQ3DnokBpmFUwq/aN1SiGVg2U59QMt3gfNDNyYdaXM90KhFrM5pyfvDsvcDDo8NKSM9gq9r6y5UyCxXEMDeJKcjNvB7OKobgvXx+FXmtoayLeA4fadJabWZeaUenIoSGqr2tuq2A0QOrRLky7Gtqx9gv+2PQqttBa9B+z/ntHZWgCR+NJsQuCjPNRSccBp4XJsqRDhXQaUNQ8m9t2ck8tmefEmOYNRehWtJtwBSoPQeXc3kgdag2RNuE1CT0QiHuMSJrF9QiG6USkAGZPL3Fd+s5MumqptS2PCWVHKDHgMyn7NTeDvoejq7jZBouRjajp0FbgUJcNhI1xQdJgMR/aKUhsHW4rQzPoEezA7WwmzIjdcaWgcQOkc8XtSpHu8hbVtsYWvUGNIGHmYECSVU8sDdIk2msOlpqD4CvZtpkdWJzXjay9K5pYBr6xoHAkYg80IgqfS1QuASiqQ/O7moNKY/cqAEBl0DOz0FmZT31+MdOGtqu1Ky+0i1CHsGjupxbvpeYm5xNzs0xtbW1XMlAvcFiLhaF8kh2uwvL9iDmMOd3bWhr6LqbrdnLCzEr/Zc0oNft3OFFk3NYBoCqBatlQxdLD0QoKSR37xsQrYh6nrcJXQsl1NiEhISFhGASMWlzrrOBMCQufOpwUaLkPIPJ4ssfimYrbTyTW0CI5I9MAC4O3CjULcuXKOixCysouoCT39c2ABA3GHRqN66Pc0hTgSSGVUZQPKAyEM8+qoTR82cwohdeWDPqT9Zjr2nPlgM6JPiYZYCXTqWciuC631FKgidiNvcIYQzvz/AysPJQRgVXuA8p0pBGR1ZjEnuaXnJ3rCrS4AlpazaIWmqqlQtT+w6bO6R707i3m0GTPGKdjZwUIisu1tUejdNc5dq7WjIwJNTXnNVqHcTxY1u0yYrQCEjnUIjJqtAiO2p1F78lpaRTdh4sxAkzsxCRr6MKJeJ/ueQbv1x8L0/m0tIehbbn8wFZzQyUa6jqE6UQEGwzlPuZSCgoC55NmU1eNh4QybrO+kwbeUIZt9l2/Fq65aK5vChY2iy7KKvaFBdodPRZWUshBCBBJf8F8pI6+Yg49h2QgmObQfVJRuKKfpDQAyJVZ201Ht1suUXt1M7nuBtDQG0TARLWFoTvmkCsywr9rwFgxCBCV/rllShhEomcb7AcCOxdYg5aHzeAm+1DUDppfQbY8ML93bwHtXICyKfErCOYS4cCq7IqPvjixrooS7coyZXJvCW8oaddibp7/fhm2MXi+fp8T3B2PIhBI3ZJHlqkopKKUeLf+3TvX5Lr0EzpyFFTX8+3K0ND3zqM1KzifATZPFVdl9zUbIMVZJCQkJCQMghRBTVJuqOsXHBEMMkCPI01DQmszbZMZa10RAABhuJZGa19PD/UUU1YqD9NKxJTYEJ0VeFk1bQGFS9A7MyMhpKmyIDWJsVo3q641njrudyf1AHPML+hky3TUl3DcCRZlAoynjJzFGu+ZRrNws+BcGQOznAy7+gqu2kZPf0vRO+nwqyf3Prqec/S8Sdu+Ipwi3DrvQUbTHhBr75Kk9h8279Cmx8/zKTSFA41bqAdyVTe73bRRaCAqhyIFLTQCqTlKJ4VJlvnfh5UOggBdsZJQce+lK6u502Q0wmu64nKa1DaNM0nm3lMljNiVXRzJpemQXk5y0SKxiNEQ2F3nNEjWXj3epvE7aRbXMzo6QuCNJIVH5JVEugoHenm9pKVaH9CAIJFCJMsBXYXUVNeg13EvFAkLf4qltXzTogBBIKabGhuBCRBsbknSTpn0CKKGi1ZkqI54YKjthygpCmZDfzhal7iJ8M3IcNVy1bVMCAu5rCvQBPTR8rBfQOgR1FTfOEGhTYq0NoLBltG4aFbhtoTg0+P6SVfIrj4InphsAnrvNpBz6/Yvr8OWFk9+gqV/jRBxlJqS7t2KgohuSRUWSqHUzbvykee2oky8b2m3Eg+raQ6ac00b2o/F3X/Lw0naIqJj3c/XCADuC66TQsQKCj85qJeA7v7uNwY1lO+NgLMnLBISEhJOAJQiuK9juFQegW4cBtO1Zp1ylqF1syk1EFtOZ3yEPdardYhUIQSrXfTMhk08B1rHPDooFKdxBGEe4nitG0O1ibMwv532nDUW5ybruL+c/Ga8KFRDJ9l63AROh3EAe4UKDJqN8dtk2HXtMIsNWe1B18hU5j25ckWYLK/Ym7CzbvFsBr1iWm5cPYZsrhtqr9KhITXqMxS/v6iOFrUor1vY961yYHbeHFLtz9DtM1TWsrnfugJl4nxS4YzZe84ZelIJLZq94wMhVybuBDB9ZF6zj2vJlEgvYqlE6ZMRrlLY7hvNMWo8yeLnqasmXkLSbi0DtoiFqsvQON1BR3Hd7FOTWUND6Ro+r34+bV23EQgpzuK6BLNXUYPXl4uMPwMeNM1gH3khSSgVnitANKDiOmHl6u6kodi76Hbfnw64AD/A5BNLRZFvG9sPhBlgcODl5Fe4U2agds+KiHwyN+mlY0+N7jVcU0CuKbFEOLDIBH4ud1LTlmagyRQZLyB7r6QrzNyAWc2jgUb3U009wmGQp44HqFaZHQJC1ifpwj4ho2uwLkG5tWEcXoJ2/WCy05rUBCAF0sumPuGxJ88nETjqy3P9hMj3UbZChAUlWWodLKPrc30hfG+mqEY4xJ5cTjiwFRSdAbBOUPQE3hkBaY/VS3DZ2I6kMGghmlwZwWL7lK4B6wkZuNIfBURQxdkZQlfh1O+UiO4BcAVG7FfM/AwiuhXAzwB4IoB7ADyXmR85rTYmJCQkxCBK6T5OAx/PzA+K7ecD+D/M/F1E9Hy7/a8HS2AGL+eAyswMPrfGwwH/+S4NIZxRtusAOjQOIPR6ii+TwUBB/IZokqOh3LR8IK03mP1M11BXGdj675cMkDBom5m/+V0L75haM7JMLnrTWDhdWhCpYcePQsZHkK5BVpspgCAQMDZ+unxFy5qDtpGcoVt/em9gFjNPT130aRPxLL9jf+++aHtQmxDXtOqJNAqz2xpnnVKga1Bp6RCVGXpSBl/G2qfre/VSGLGt5uX6o0ZwHZMCubgiatLvE0wuM/9mSOF8IehRhqevmsW9mvfdPB+rLVTNe2uKdF5Itt2xh1MHLWXKNFqF00K4XIAXh811KzyhODZyO6gs8m7cBtLiR9cCPhPAs+zvlwL4TawQFswMLkvQVHKUsIS74G17vIqagjoGm8i7CEDbo8qpx9H1LG0dVJk06PJaqX6rZsAI2mn556YdMS+kUNkvYFlzyw7QuFM2QXlDThxOUOT+69KiXY0rJCAHbxsoRuQDAaUnjmtPk8I6DMoD183Awl2DSTwgDwiIVUKiY3+LpmoF2q0QEPK3GKRYRwOpEx4A1MJkGdDZBKjLIC8ZS7uEHhhYe/qy73duUpHlYZtbbsLRIOwnA7GnVkgnoY5ciN2kKbbDsPbnOs8oGbUtvx/Slaee+HAf7Nxq10wxHnhNqXqUy+1aICRh4UBEnzOijDkz/8oR2sAAfp3M6PGjzHw3gPdj5nsBgJnvJaJH9bTveQCeBwBPeHTnKQkJCQnHBGpcc28ArNIsfgzALwIYmIfi4wAcRVh8LDO/zwqEVxHR28ZeaAXL3QDwjKd8IFNRCI1CaBOWokGWhR5RVRnSOXI235rZigAfiBlXZ8PETFfuVzlsDo5mn25mbMi0oBbEqn29ya2bdjrjdK0k7WONjH67ucxrHdyc26ThtjRJ3Rgn5axXLiDltA73HDPxXJQ1tvs66saTZWJjDEhqE9E9BdRTTDsNaYA9s/5eI3efFjGG0uqjnMS2+y2Ns5RlnjJRdQm9e3NzrtM0gIaakukvXDlWU+2C95Jj4VVFltoR8STy3ObiAXo21hCEptFK3S9TqkhtkHWQrZfqZePlpGuwrsELk15FL+ed97c2dA0o25/yYsXJI0HmPd4oWCUsfpWZv3zoBCL6iaM0gJnfZ//fT0QvB/BMAPcR0Z1Wq7gTwP0jygGXpQmoUplfx4LzafBBkaA0pLpLSgH5NLQvOEHiOqzqGGwiIeL3SdjjpGpQphtXx+AGNJhzENlrsyz0eok+RC9IqqXh90vD687yGdgFe3lBYOsXwkIRGYHQbgnIcsMyitjTT25/vWwuyJp8W8HgmU0AUlD2XOmFogSt4a6jaEAKXVDDAbwlIOTvLg+oIVopKmfYFbctCAK6o4cz57oO2sCAH/S5XEKpzOQvAswkxnm35ZNmsuPaNyK6G4ARBu5clQcTERA1848eO1pQZzxxgvhGbB1K0oVKhe2WNKquoJaH4INLZnMxb2wN9vqNqCel+s9XKrRnbAFEhOwG8oYa7CHM/EWrChhzTh+IaI+IzrvfAD4ZwJsBvALAl9jTvgRGu0lISEi4pkCZGvV3FjBaLBLRx8C4svprmPn/O2L97wfg5XYWlQP4KWb+NSL6YwA/S0RfAeBdAP7R6qLYzvKVUTNl4JqgUwKvjLr0nhYMQO2gScHA2msUXC6NuukprizUJnRjPAtmkEpB7d3UHFvOgVw3PGc+BcRsm+qyWdmvhvCRD6kIzgrvMw7WQLVosqJmE2gX7JXPACJk1ltGpoMgXZo8PGMMwKQCn3gq5+GsWbRRGsMN3bAAlTbzarUMg816DK6oK09xxO1pGb/jY8Cg55KvbsgI3qE9NIdiQ3hbkwhiAUZoHVAZuCoDatNrjvZ5BrP0WryL+FkEbRXOFVmjAXRpqmNTYJDQJMC6SRsOhN8WAChhtFcqiJ3Qlx+GPrhsiimjLLAdaTs6Ee8fisEQx1r1bQpK3lAtENH/APBkAG+AD4MEAziSsGDmdwL4ux37HwLwCWsWBtY1CAVY6yiFtGmyPtxHfflhwHKW2YXbQEXzERr3Rjt46xpcWlW4Ks1CfFZ4kExzXtiB21FNWDZBbdMdcNUIJMBQUUFwnxtMnNqtbEdWGUj8ZtsGAMAspNWYVEMFsW6WpHQDt4ggFzcL6WUUCw3O8kZYcfM8UVdNsjeYQY3qKKeWY0VYgxZXPUUG1k076xLey8s+C5KDdExL9VFMqwLk5P3K/66e4JEMCIK+a+IBqkc4DFEgFJcbeMKFFBOr3PPk7N5Lj3AjlTXX6zhgT3KSA8Ii8oaSAXNGWDSr2nHk2kp50Ux+WIPnxvur3r8Cfbjv+3xnric5MXN9NDYmDwmSIcqp75o1QV1tOsMYq1k8A8AHM8dEb0JCQsINiqRZdOLNAB4N4N5jbMvR4AKQtIa+8ojXAvT+lYZCme2ZlAGH+/aSDGrPUDae1nGzpMXcG9nIzm5I5JmRM8UgB01eAHYG5curGg0Buga75cu09m3x5zoNpSiAyczfG6Z7/hhPdr0Rm5Xxn2fvOVSKYLa6PWN0s9DYfx7NTJxdvAeJc50HTF2a2aR47hLEuqliecl4vTgaME5HIa+L26irUHvoM3DLe1jH+Ny17a7pojPGZDpFhwYRp6DoClyL1n2XWppzJpD5o+KVE6UnUYiGBiLn+NFUHv5ubYs3omv/zrkqm2+kXAbPl6tlcL+UF6Cp6cO8mKO+ctH8Xs7Nd9jzTIPZuvgmWGob/ni/AZyUAk137DVCi9+WNpDiLBoQ0S/B0E3nAfwZEf0RAE+yM/NnHG/z1oHNDVWV0FcegT6wieci/pMmM7+veuC9UIc3AwDy2x/tB3N/nfPMgP2wncud1o0AsJ3Ze7IAoB37EUdrZ7BdBawZMPaD42ZNDSvYikIIhx3oyTnwdM9WEAbzGc8aMZhb+wXVNTgvwgApyVt3JLDz96ArsE/vXHkBRNXc2B7yyENHHFfSW0dG42oFJZ4JqzwUXrHg6vDAYfFOzP2ssBn4YwNeNetEBY8ot/f8WLDY/4Qo24B4FlRWpq8V0+ZYH2Uk3XO1DgdF6eLpBIhwISWgoYyi8nl+FbpsvJO8Lc8KB0/XdngBOpqXF4fiumGbAcf0UkxDRfaeXigFmpn1QyifbE9IWNxo3lCr7vR7TqQVCQkJCdchkmZhwcyvAQAiejEzB+k2iOjFAF5zjG1bD6yhl3PoSw81cRHxKVVpjN92RkVKga1Hht7fM7MQO+uJ6QJvyHb7rHpNeQHOZ2YG7465mbXNkunK4uXcZ4SNQZ4ucHTWxGsrenoePD0HtqmVjYfR3P8OAth0FawGyCgaDySZFpoZyKq2L7x9lqDMrIcMqzW4WAnn/SLqY8B7PKn5lZBqYu0DzCgrwIdG4/MxMR1BTVzXZtY9RCF1oHNGv4af/kZ++GtqHb31LuZixh95lJHyyVORFw1NlWVBzA6XgiJymrFzpkBEeykVtj2fgBDSpQCgr1yEnu83WnYtPP+WcwSLEQ0+i3o4Y6yF7w8d2kNL4wCGn7HKgrGAIqr5yEg2i058Etq5mT6lY9/poa6gLz0EfeVioG5SljWd1LnMFe0IzvrSQ1C6bmwIwp5B05n5mAqRB1/yy1njdULV0g+k7uN1ajcf7rfdbi3Y1kOqEUK1dYHVu7cYm8zStE0dXjKDMgBeHhp6zQmrfNrYM/IJuNgB26U8mTXICZzlfugGKXlyAKh1445blWEgWNBwkxKebJ4jnu8Ld8kscCumWdZE5s4jCk5SgJtizdxB28CYAdAjTmznaRVtKB3nNFcUxr4GeIHrbQZo6CMmZd6V8ypaNna2pk4RJCcnJEoF74mAJrV9uQ+eG+Gv9y8HdoqgP5ehjWL0+xugjrjLruBcbrMVfSQql5fzYBLnBd4Wg/OSN5QFEX01gK8B8GQiepM4dB7A7x5nwxISEhKuZRARVEr34fFTAH4VwH+ESRPucIWZHz62Vm0A1tyonHFaBTlb0boJyokNjko1sRSTWaB6U44mXw8g4ipg4wfCGAHAGKn1oZjp5ROjBTg/+boO1WTR8TibQO/dBgC4VOfYLRSmlnpSh5fAiwPfDprMvPFbT897TcKlceDCrvsMQjYx96BIAcv9xvgcPkxDWS1twOLhFWg7u1M7xivLe0TV5nxtPV1avu/RM3YaBR/uG0rwtGZmW/K1B7C5RqTF+5bUqa6bhXryIjQIVyVAzoPPard1Qwt5Y3MX3D1nWaBpqMkMXAKO6+LlHPWlh5oyF4fNNyQcNLq8mUa9T92ktB86J243x0b7uO6OMpoxoQar5vdWQICaHN3ATUSPh4lZezQADeBuZv5P19q6PqtsFpeI6AqAD2Xmvz6hNm0G1t7dlfKJsTG4Q+7HALfKZQl99WKgmmvnYqtrUF9CM+epYb0u1N5NzTE7kLsPg1VmaALpjSLV+7IJcDIeVoYyurpfoWag2LnZF0279sMmZbyldm4BABxUTRryaUbIyQgJAFjWGkTmlc92bzG5fCy1FXgfaRN57QQSi9w9uipBO3u9g61ezoOBgPLCe5HxYm5cmQHvOtksB7sFobGWAGgL/iNjDYqEICgs51Xk6B2UPg+Po1H8REJlIkCSDI0oaKEmj1lUn7BRcAVLq9o6LFXo61/MvS2PF3PoqgTXI+0v9rxVXH6f22zQXgdHY3ZRleI+O4PAhAeWK3OVN9Z4bC3rbAXgXzHz6236o9cR0asAfCnWXdfnGLHyTtk48L+RiJ5wAu1JSEhIuC5ANs7iqLmhmPleZn69/X0FwFsBPBZmXZ+X2tNeCuCzju9uVmOsDnUngLfYOAtvmby24izI52KS/uOG6rAzqsPQqMpSHdY1eFE3QUTLeTMDuVwHBknWGspqElyV4OXc0FYAslseheyWO0yL8sLM9tyMzWoVjVeKAp2/2Zxr4zjUjqGTWCmoQ6Nx3jy9GQelxsPWH3+68yjMLAVRasZhpXH1qpltVpr9gkLnJgq7hULhM0Y0irrOCtBktwngq+aNllSXdoUy4RcvZlA+V5bbruvGsFqVYDdrzwsQCn9ufflhn/LBecd4Sg4rMGr2H80YR2sM5fHTYVFAGQPe28zRQo2WqVELRwvKMt+/aDJrvPl0FfTTIGDOzaTjgDYLAryjB1eleW+WwtIHV7xWXc9HBs+5fZmjxoY1h/h4PKDK477MDtpyaO0EAL4Pm+vr4P+RsZ431O1E9FqxfbddYiEqkp4I4GkA/hAj1/U5KYwVFv/+WFuxBZBSZqC1dgE/0ItF3rs6CUsVF4IakF4unh6wdUmqyarw2nK8LpkhAKjztxhvLFdXaReQtwJJSXfbnT3g/G3Q1l2Wqjnyh98FALhp5xJmt96F91017Xj4kDG1y9hlZFa1rP2Sr4Drv6VmLCrt17hQ1AgMZmMXcZ5TJtLacdaHYRCVblKZMwAs64Cu4+U85NzdOgl5Ab2cg13k7sHl5uMd4L27IWmjkUnvxhofVbZVD5m+OoBwcPMCsrSTCDd4l0svVN3zUbvGTqEAaPd8M+NdFvDycl0IlaGLbiOlwDoz7rquLXbSAwD6yiOo50ZwuP9yUJSDtwQp1SskVr1neXxIcHOtOwfoFkUV03BFAThPwy3aq9aYZDzIzM8YLIvoHICfA/D1zHyZ4pTzp4xRwoKZX0NE7wfgI+yuP2LmlWtMJCQkJJxZEG1N8BBRASMofpKZf97uXntdn+PE2KyzzwXwEpi1sAnA/0NE38TM/+sY27YeiEKtQqZ+kN5RHbMhs7/2OaBa6DCs9c2UTJyHNPBOgqA8GXegheeQUhnoPHzAlTq8jOU9zaKB0w9c4qab/zYA4G/KCgtr49wtCDs5+TxBRnuwzbbX1oLjyYKprcgSWpVNuvRqGfrsaw1W4ayNxYxVqvoyTYm+9BD0wZVAm9DCuCiNrmPRaeTswUqtxaczOR6tItBsXFoYd0ycxxWMt530jnKeeNYTSYv4CW/8hqE2fZp9QSWZflx2a2GFyaEWZH3Vtc/dVM+XqA4X/hgp5eukrK09SE1jXUO4vF4ek/vcs2hRUJmKNJK2F54Pnp3MgqDa7YCaFEBHKcWoEP8NwFuZ+fvEoVfArOfzXbgG1vUZS0O9EMBHOG2CiO4A8GoA146wYDaDXFUGSQAlyNoL5IDl1fy8CFXXIcpKDPjxsXq+BC3e58tQF24LvJ94WTcUgxRq+5eR6dp7U1WP3O+prYN778O5aolbPtwcm8/uwGFl2p0pwk6hvBDQcvkGu+2OKYJfOlXBrO2hFlfNjv1HGrrMCVjpOixvVkYGl2VjhwDAczT3F3nSDA0QYzGetoK3VfUfPwYfeRUN+DJAFGjet8pA4utjGfPocpjBCAuuls3yojIdfpYFNjGZzyxOdgkIqtW2Sz7L8vJl6KVpRF2WvV5NXcJgjIDQK9aaGIpXiG0Wff2olaQxHsi3TTVub1nVjwXwxQD+lIjeYPd9C4yQWHNdn+PDWGGhItrpIYzwpEpISEg4u6CtuF0z8++g31b/CUeuYEsYKyx+jYheCeBldvvzAfzK8TRpM7DWxoPDqtReja81lMwMqTKoiaVJ4tgJGVAmZ4hWS3GeSpKWMqp3M2NhrVHuu7wN96EAoG4xTgzkPJ3czHu+H1A4en7gDZn64ArqQxPnUM+XuPqOd+CmW14PAHj0+38kHp6auIplzVjW7LWHTFHQ6wpl9jm4X1QvQfMr4Iv3mfquXmw8bqzBN7ivYPU3kR+oMpqFDHSsF4bC0GU16PUylrLYFKu1l2rkeSvqiWICuugQSUEBNs7CaRNKg/K2duF/lhV0ade5nuz7bMAcO3No2Q/rlhbW0DmhcwHXps+696HLKrgnUmpYM+jBOlpgfGbXG3G0U5dROaaqdFUaBsBRT7GxfxsgbEVYXC8Ya+D+JiL6XBh1iWDcvl5+rC1bF7pGvW8olZhX9R0pSiSmJjMvCDzF5C6KuHG5Ol6g4kvVH9Yry6qm1f4clD2A3F7rck0FwU9yzYxqCX3ZUE/V5Uuo9ue+zHq+xOGf/QkAYKcqcfuTnwYAOLz5Cbg0r/3yhZmgnXJFmGQKWW35Z7dOBQCqFlCLq2blQFhh4AWlTcFeRdHn7nmKY7ycQx/u+8FMl5WnM3REL3Tx3ccqMMqRgqBcPyiwS+jFXkIUU1ACHO3jStBm0uZVTMC1biiiw4OGh++y34hJjF5W3psJaKKNnT3C78+y4Dxd69BG0WOLGOsmuy6kcFKBi3bbE6rv2bvnpSaWZo7sZVtp5/aC8q4LjI5VZ+afg7HWJyQkJCQ4p5obBGO9oT4HwIsBPApmQkQAmJlvGrzwBKFrjeXlfagsg5rk4azOzUZio7RNDQK0KSmCDOIxC76wpYUoUz4orysTppxxVftzkDILDPJ83wb/ifgFMYus96/4Gd7y8n5r1jR/6BIAoHrj67D7iDEh7X3wMzG99Qmod2811wnXpxkqZBffC1raDKKzm3wOKVrsB4ZUqUVwtWytgkbyuKDPuFqiOpgH3jP+moi+kHMwT3kckydSF1YlfRtLR7W9cJr3PVRGl3Yh98u0NPGMtaE2gWzPaqPu3QnPPzd7Li8fBIZqAMA81ChcW8v5EvWyCvqtu4pU2znBt/+IGkWfdhZTTUPG9j7thjIV0KDV8nAtWmz0PSTNooXvBvDpzPzW42zMkcDsuV1V18gmxhNCFXnjjaPCwVmXFbIezjH2JNHLKvhgJzfppvzYHVd04HpZoraDPOx/1zZ5rqujmjd8v7btzuxrctvzhy5jcdG8ir0H7sf0sU9AcecTTbsuPMrnE1IHF1Hd924fjZvf8VioW+9s6mbtc1phqcJEjDq0UbCwUej9y552qudLVPOFV/tZRzYigVq6MA8MNOsOQmMHedfm3mssS9E1AKyiQ9ZFKzhPRLN3BQmW+26CoVHcZPOQ+RxdLgC1DDyaJLUENP3Hwa3yVtvn0uW1xpkO+rcsQx3hGUivpnWfpaT5hoIAudaobJ+uyyZKf2uZYmk7Bu7rBWPf0n3HISiI6PFE9BtE9FYiegsRfZ3d/21E9F4ieoP9e862605ISEg4GqywGPN3BjBWs3gtEf0MgF9AuAb3z/deMQ592RYB4PuZefSyrmw1CwCt5F3xTNep6kEKA6XCmWc0u6znhm4BgHpZ+dlNPptCFXlvCgQAgfFXntM1K/LnimN1WQWzONban3flXfdh/96HsHPbOwEAk5vPNfmlbAp036ZqiUx4dsn07SagK/KqEdqEDyTcv4zyyoHXgOr5ErrWwSzRURfxPXY+m20YR1spoVbPgagen9eIMhVoRaSUD+QbS211GojdDxew51PpCy0DxjA9sdrE8vKBpyN3Z7smLkeFTgSuPvlehmikLiNxQLP1eLDVK/I7xYj78Fj6rg9d/cm1r16WwTla0mzboqO2F2dxXWCssLgJwAGATxb7GMCRhIVNkuUSZV0hIpdtcSM0tE0z2Oe7O8LLRzWRre4aJziUioRHE7xnPry6KX+SexWfa41ib6fXK0/ypmQDody11WFIE8RCznsVlVXrYwrOKyvs/43xotr/m4eQzYwdZnbzef8bsNHl1kbD8wMTtOQEQpRjCBBeW+US2q6ctnjkClhrLC8f+PuXH77k82N6Q+I4vaC4Z+AJnmEHndKXn2hoMNNoq+dDg1JvXqXovDitfj4z6ernD13G4hHj9Te9+SpoutMkztzZg5Iut2JSE4NrDbJLsrqJkrLDgXw3GlWnYB2LoJ/aco9CX63qN24Sw1ZQuj6oegTe0UCj85SdBYx1nf2yoeNE9AJm/o9HaUiUbfFjAXwtEf1TAK+F0T5ai34Q0fMAPA8AHnvzuaNUn5CQkLAeiEKX+jOOoy/zZPCPYFbT2wgd2RZ/BMCLYCZYLwLwvQC+PL7Opvi9GwD+7uMe5SdjbrbbQlVCV6XP9ipTUZBS4WxSzNDq+dLMvMTspF42dBFlClnXut5lGRiquTZahZvtS40hpr3kPdRlFSVdVd5n3mlEnlLQjU9+NV+0Yh1y1ayAJhdfMqusheqRiwGpDuaBcb+eLwPNKtaG5JNvBYYdo0bRNWvv8qqR5wINnRJrHmNnwPKO+q7oMsb2eQDJtCBxu1XRaLVX3/sALuyd8555MiBQl9VgriY383b1u2vi+tx2n2a16vluasR2kP13jFYhqWT3jQLmHbs2HEWzaSFpFmtj41y6XdkWmfk+cfzHAPzyyoI44scdd7lYeJuFtznABcIpP3C7qGVp95CDruOAXdnuvCyb+PMAICsK7w6qlxXqZRV8hHVZedU4K/KAvgnaPV/6AToecGsAqiz8dZmwyUgqy3mH+ahW6ca6nEPlRbMCnqDL3CCzvGzXNBCBdu54LT7KLpfY+LeD3hJfrLqieDsGpqMMVmOEQNc1Rx0+ugSsD1Kb5L4/VftzLO5/EDsuZb6Y5Tp7hXdR1jp4Zlpr6GU4OXD8e/CsouvC87tdiF2bj+IxJj2e+tyrVZYFx5zHmGsbEApANw700ZRrY4tZZ68HbEtYrFy7pgt92RZdWl67+dkA3nz0JiYkJCRsEzeW6+xpaxZ92Ra/kIieCiOE7gHwlasKcoZjF1chaRM/uxJxDIAxGkoDsJwlSe1BL402IL1H3Cylni/BtcZk4mIhGu3EXeMN1dGMUdbtqKQxtIGuNdSy0YDqTEFZGiyfTaDRGMalpxZlyj+XrCxA2TLwHvHUkjZalDSwyxmqXpbBLFSLDLVxm48LfSWrAd/7IXRlMw0NvaKONdp2ZC3D9sMuWrWaL3H40CXM7rTrqM/2fCqQfDaxgXluZl2jElqzimhT6XihJjkya1B3AXryPvoopr68Teveb1d6cl+31I6EVlEdzG16HXO8YQzaZUmq6kggJBpqA/zPTS4ayLa4fpJCNoMzZcq4mlpPj3L/EJOb9sJTbScv9w9958r3ZgHVosvKu99V1mYRXF+GZblrA+HjPIUcVwrTxmJv1m5+xwDnPGC01qjnC1TCTqBdAJkte3Le/K/LCkp4r+hl5dtUAVDCjbcdENiUDzSrpBmhFyZLbLd/S2tMYHVULGVZb66iWCA7CqU1sK0pULJJExiphT1j1f34o5Z7H6LFOl2Ll5WZgCzb9gTKTM6wxf0PAACmj26WX53ddgFLISycizMAFHuzzvflqRvxHTih4o7VZdUSHLHtBWiEr3zevnz7/LqC62Ih3Yra7rGBkVKo5guoSdF5nbyvekvCgkA3lOvsKLFIRN9NRDcRUUFE/4eIHiSiL3LHmfk7j6+JCQkJCdcgXG6oMX9nAGM1i09m5m8mos8G8B4Y76ffAPATx9ayNcFgT9VIjxwSwThZUUAVOZbzJkbAoV6WPq+UOyaD6aoodYK/br5Aphs6K9Ys4llynk16PUiAZubEtfaz5Jie4ppFzAcCTUXm9enydnH0mXsuSsz2pPYkKbuYnpMzyOPIt9MZnxA8M5ENOErDYfaFhlyHmKLqCzaLDbdA6FkDdFNNQ+kn+q5xdfcGvtn8Tl2BnFxrlMtD7N/7kG93cfPNAIBsZxeTm3axuHiltz0xJF3Z1GGftWCt3DNVkfYgsSklNfQs+s539bHWqJxTiPdSbL5nScltCyk3VBuuqzwHwMuY+eFrbTFxkxuqRAUgEx0sm4XeSnop+PeywtLxuPtzTG7aDR6ItFHoZdmsqlcUIS0j+H05yOQ7k9YgI6PJuwRFdwBblKcqcqXs/LjFtrNTUJb5Z+O8pnRHZ3eCom+A6ko1vg0Mec/0La3Jdd2iAvqontgjaEhw9KXCdufKOnywWW/r2+V2RYTHQtn9Z619P5J9xCTLK3H4gAlBqudL7N5p3nU+m4JrjdltFwAYz6la2CxanlCqJzmi1qAsCwIs/eTDPs9B2qiDXnLPTFJ5Y1ycu457O0StUR0ugz7r7Hnx/W2t795guaHGCotfIqK3ATgE8DV2WdX5imsSEhISzjYoaRYBmPn5RPRiAJeZuSaiAwCfebxNWxPsKJgSXNfeOMxaewqJopXbuG6OZUWOaiI8h6KZlvQIUkXhDWn1fAEWHlDmeKP6yrgHSW113oL0vhEzSGckb2Z7jHzWraEYKiQ0RjcLwYjXXYYxGlLjqaOALhljktn07/LYENbxtd80JsLdb5eG0aVdAO04jS76o8srZ4hqcnsDiinSiMbMauMsr7LftLQ8cR+Li1c8dagmBZRIYaMmeRDIuVi289OEaf2z4P/Y2Xhs7B6b/6lfc+zQtIU2Ic8zfb3223Hut60G45lSk7CIQUS7AP4ZgCfApNd4DIAPxJhguROG9+SJ+HegI/CtrKBLl/+pNi6DjtPPVKP61zpIHpjVNXJrm1DWC0N6RXm7h25/6KsGwnAwsIFYReHbCQBZkfnBPf4AYpWeRLR3WE8dlCnr5rpu2UnkOZz1u/gO3c9RgrQ2xRAl5RB7TLnzu9o+ZvDviwgf08542/Uh378CId5eR8V5sJX7c2TCBTafTZpJjNaY3LTny+mzx61qq/M+Ux0CYaivu2Nd6c7H2MBkcKz03JOTGtTGntfVhq31Q4JZffIGwdin9t8BLAF8jN1+D4D/cCwtSkhISLguQMZuMebvDGCsWHwyM38+EX0hADDzIV17Fu6VM3cXJOdmZAG1NMkDWsodB9rBdGpSIO+IlQDMjE6VwgNjvuyd3bpz5G/dYYTjujapSYTGMqRSx4bKPk8iAC1Dp7s+oMTi39GqaieJIS+bLmN3c2wgc2xHSotVWkZnrMmAoXYIfY4Esn1dGueQwVZZzdjnMCtLr2UACCiqTeE1hMg7CujWxLqM3W57bKr6+Bm58suDQ7BYJVJZ+jcM4hvniLAWkjdUC0si2gFMWg8iejLEuhbXEryKrpzLn/AW0hpaLDUpPZUUmsG9VZbz3HABfLMJZjef9+fpZYXlFbt0qdYmmTu68xd1tdX9jimGvhX4ZDr1la6FA4n85KAno7L9ffVx87VGtsFAM4aGG3PeUQSGwxAtBaDXa2oVNhGekmLqKqfl4dZzXgwlBuHqcOlzJ5FSKPZm3q43hKzIe9vltmPvqK5jsf0C2CxAUgoISUOZugi5tSU6OnkbadF72wKAk82ihX8H4NcAPJ6IfhImTceXHlejEhISEq55UDJwt8DMryKi1wP4KJj0HF/HzA8ea8vWBo3y/yelwnQGYrYOhIbxeCbqDITZbOJpKMoUJvMFyoPGkziLYin6goyGqAYJV0Zs9HT1x2UF9xvVHdcvjd+e3pCGwh5sSj+N9Xja1HvGXNuVyqI7FiM8pzuYT2LT3FOrym3P2IfTq/QhfmbBVod3V1997rku5wvkezteYyCxpr3rs7HBu1XvlhDTblyLxY2K3BjxhccXYDwg5TVbRxIWIax94lMAPImZv52InkBEz2TmPzre5q0PZV07u6gIpRQwaUJRKXJjlWq0+RBCbti5y+azCbJpo8LnsymmdvEl6Q3ltuseb5NYvY8/2ECYCAG3TtRor5un2x95inVd0xlRLSPdo/PHqPxjhMYYT6qhRHNhWZsJkBhHGXL6JwPtto0VEKue39Bx1//y2cSkzrdutyYflT22MzG2PGvbil3E4wh6SUkB2xUaTjCF303jNi2zJzjX4Tpz63uUxyAwKHlDdeCHAXw0gC+021cA/OdjaVFCQkLC9QJS4/7OAMaKxY9k5qcT0Z8AADM/QkTXVnYsoibvUY+BU01y0IhZM2C0EHfEzbRym7xWFblZvxoAtEY2m/j8TLoQq3NFqaVjCiLWJmR7uugqaazzCyihiQvpw8aU0UDacelFpqO1xNfRNFbRQX3njdU0WmUPeIaNQTu2YTjdwzaz8a6L2BtI1sV1s/iRf5cibXlehIZidxc8X/q+lwPIJ0VnX21pGEeg78bQtS0WwDqBZGKI02h7/h0JZ8gtdgzGCouSiDI03lB34Gga+bGh6+OVXkQyl1KFxnWUsixwJzQqdTdNwVoH6zeYVe2MzSJOc97lDjvER/fZNGLkUcLCLsT7ZT6eoXO18DjpOi6j0p1Lb5wDK663C2PaEuMoLrt9Lq9HwboJ8zatf9sJ60xeKcv3w7TLCYvJ+b0mMeYkR7k/D9aBiaP+3USp6vAkjIUG0B0EOYRVdK0sy0+sZrqVSp3qY5jhnxGtYQzGCosfBPByAI8iou8A8HkAvvXYWpWQkJBwHSC5zgoQkQLwVwC+GcAnwHhDfRYzv/WY27Ym2Kcp5rpZyU5nys8wvKeEnWHEXktat2cjgPF+IrEAULk/R3X1qim/rFDuH6KymoVZUCkLyog1ijF+8n2BVr13P3aWZs8r9ma9hve4TCXWCnfPV2oSQ1rFmLbEOMp9bouOWyun1VhjtOqmhY6CTWgwl/dMpq43/dQZixV2bLZaXdco9+fNmvOziY/PiL0PXeob2bYuR4Z11mEfCkaMj8t66vkCiOJImjxX20pRTikoT4KZNRF9LzN/NIC3nUCbNgSF0ZpOxRacaucA7LyfOkp0x9RsijrLfCdbXLzqz8lnE1TzpU8DrsVKYu5j6aOdugKcxmKsC6qutY/8ziL321VwVJeso5ovAzonFhTbWFpzE+8qX/8ID6o+bBLRvOm9bnOFtSFPttBbqXuQjO87yCe2jFbOi7IX1PNlSD85m0G8bOuAC/kQVg3sfbSkCbBd+HZUYlXMTSc37coJSN5QLfw6EX3utZfiIyEhIeEUkbyhWviXAPYAVEQ0h6GimJlvOraWrQ0OcyJJddRCFXk79XZkIHMlSGM3ZWamJdfwddqF3jV0jpvBKTSzObfIzLa0ic67Hri+NUMvcnDWr2VJyNmlKgq/AllvO+Q9jox9WIV1UzSs0gzGag7rtHeVhjCmzlVpYdaCiBOS6HLY6IPpt6avL6/sB1mVXVlAewVJGcugYXJRjXXYiLEuVdSb98tmVs5mE98vi568bpsg2SwiMPP51Wddu4gjn912GK1dBxyrBnz+I+dB5TxE9MWr3otq2aH2n1aSPYdYSNQi8M6p5XW9aJ0brPI3m/hzpfdT7H3oaCdp05BYZ+AdKxyGBuDBYyvaMmYwXTX4Dw38fdeuFHJqvBeX63NZ3AeVGogeNwNzsbdj/8+87WF5+SAIxKuF4HDruHRRf3G2gHUxZF/YmPrr6aMb4zpL90FE/wjArzHzFSL6VgBPB/AfmPn1Y64fdadE9PSOvycT0bERdkT0bCJ6OxG9g4ief1z1JCQkJGyMLaUoP6Hx7t9YQfF/AfiHAF4K4EfGXjx2sP9hGCn0p3b7QwG8EcBtRPRVzPzrazR4JWxMx38G8Ekwa2f8MRG9gpn/bOCqlakehtIFxL+DWILMrCjnaJmsyD0NVc2XqG3+HFmHK29lkNgJaiFBypBMGWN1tFASYFI8tFbci9ons3rK7KYSq2aAR9EkemfpgzP7fs1haHa/iiYaqzHEbdtU03DoDE5ToaYo+5fUNrS4nrIMWZEHhmsXN+RXmRTXOq0UMFqrfI+tANIVgbKrMBQ3tQr9+dK25VywnXQfm413G8F90J8K4EeY+ReJ6NvGXjz2Tu8B8BXM/BYAIKIPBvBNAF4E4OcBbFVYAHgmgHcw8zttfT8Ns4xr78MjosDzSUIKC6D5+GvB48bqrlJN4sB8NrECw9IySmFy064vs54vvG0k9gI5DbiPNx7Ada19BDtlyuQE6vEMkV5Nuiz9oCGDAd15Gt2q/ZBn1BhBse5g2h6Mty8cjiK4WsIjTvq3rjdWx3uOV6rrGzAVGvsc11EWgmXlXcHj1RQl2E4U5H0062e0c45t+j6Cdq9DQQ1Qb1vDdmiotce7DfFeIvpRAJ8I4MVENMUa6bvGnvhBTlAAgJV4T3M3dwx4LIB3i+332H0BiOh5RPRaInrtQ4fz+HBCQkLCsYGJRv8BuN2NVfbveaKoUePdFvBcAK8E8GxmvgjgVphJ/yiM1SzeTkQ/AuCn7fbnA/hzK5m2nHAFgPG2isGtHcx3A7gbAD7szttZ5nPqgtY6CMQjrQJVPGiAnXkDJqusjjUPl/Jgkndk31w/6GdMevV10UcPyTq74GalbnYYLBJls/rWHb71gKGnfDqV+bK3DUOpR2Q7utrZl021vVJdfx6noVxT6yyE1PXe4lXiAAQrxfVdt42330u9REZyEsGqWpwDmHY75w8tVkVcVQeLNDH1ska8cp1MH77NOJMhbNXTLAYD3BqVevEgMz+j59io8e6oYOYDAD9PRI8ioifY3aNj58YKiy8F8DUAvh7mxn4HwDfCCIqPH1vZGngPgMeL7ccBeN+YC1WH6jmYflukOJbnKmGjcB+WPwb4wTKmGoYExZjo4+MQGBJu4M93jCB01IPjqIGOACqt/XmqKNqDcLR0pduanN8NVmOr5ougnrER3GMG6KNiaIW9eijQb0XOoz6bFOm2PaOfY99swBsqL468dtkPANNH1pnw9GUbqMvav+NskkHZ367vxd9cV5nxsSFBfvJg6DWkxQA2Hu/WARF9BoDvBfAYAPcDeAKMsHjKmOvHus4eEtEPA/hlZn57dPhq1zVHxB8DeH8iugvAewF8AYB/fAz1JCQkJGyMLU3/T2q8exHMAnavZuanEdHHo1l2YiXGLn70GQBeAmAC4C4ieiqAb2fmz1i/vavBzBURfS0Mv5YB+HFpM+ltp50xxdqFXNCoLkO1Op5NypmKD8JTjUbRXW/WORPbxJ97THruo0CJmW5vIFOtoQ+XncGEXVqF1ATqsvK0RoaGtoL7bZ93PV+0jOuynKGUH0dJdT1m9T2g36sq1kQl6q7ZcEc/dG3xdQ1QWa5dx5F1VraFa+37+6qUKX3aUvzOsiLz+7hmaJjy3VvPfNBn7hcWi7GqLTIFel87jwsMoNZHFxebjncboGTmh4hIEZFi5t8gohePvXidNbifCeA3AYCZ30BET1y7qWuAmX8FwK+se52kihxch8qKHOWyiUTOirw/4lXwtrLcGLpsq+x9QmJT+mTbNETfsSHhJp9Fl5eVgxzUyv051KQI7EST88aLDOd3sbxygOXl/ebagdw9q5Z5BYxAWfV8x+aOGhq8gWHBErpgN0Kl7hNWYoVGV1dLmAzUN6Z/cB3Zc2K3WuHlZtzGbbujydWYyH8AyGamzQr9tokmAFZ+Z8aNtx75ncRC4qjuyOuAAWxBVpiyNhzv1sRFIjoH4LcA/CQR3Y81bM5jhUXFzJdSaqiEhISEBlu3Qh8v3gjgAMA3APgnAC4AODf24rHC4s1E9I8BZET0/gD+BYDfW7OhxwuiXgpHzj5kLEYQeBdRCNLgyJk2Xk9illh1rGg3qplboFM2qasPPj3EJG8FVHVpGnW0bvkQPaIyhcUjVzyFlO/toMjczHOCyfldr1nI+t21Y7QJiXXOH6OFjMXYDMB9oEwF678P0WWbaBXxeU7LcO/XHXPpPgCgckGXUZyF1Hrcbx29N4/oHbpjQ/Em/pss++9vDOW0KiByK+DtaRYnhI9nZg3jAPdSACCiN429eKyw+OcAXghgAeBlMNzai9Zr5/GC7LKq3UF5zb58Ngk8cjyfXrQfhRQcocdP4+URD55rtfk4VOMNB0ApOBXCSPRa5P/pa3OXG6yOhE89X7Sun95i0o7FiQqla+26QmMMVpW5yiPLoUv4x/YGf66Sa0aEbrxOYADhEqTr5MFa1Z/iQV726WySe/fwYm/HB2HWEc06eu2O2OMqCt6T7TVJK43Noist/hiqSZbV157jAG/HG+pYQURfDePN+uRIOJwH8LtjyxnrDXUAIyxeuE4jExISEs4qGNfo2tJt/BSAXwXwHwHIvFNXmPnhsYUMCgsi+iUM0HLH5Q21EWg92sVfJlTjeDYpvXgkBSXTnMt91wI2zZsTzAKtluUpqiKHVt33NxRY1xXjwbWhncrL+2YFwmhVQcCmojgiJXVUjAka3Bb8fUdG7nWvl4jbrld4gpX7c//eZ7ddQG3pp2p/jlrUYfJKlf63DO7rakvXM4y1ADUpWmlkhp5Fp2ahVlNcx4ET7pYbgZkvAbiENdxku7BKs/ge+/9zADwawE/Y7S+EyRd1TUFlyvOkYeduXAJlZLZZIjLspFJ4+H0tbxcdqfR177nXAroEWZ/7I4BW2ukYm0aGZ1k4AOiyAvmV1fIgn1BcwtDAvW1BcpS8Vb3nb7yiXpteiTl+h3UEd/zEdFlh/tBlAMC5u56A6c3mXcwfuoxK0IfSs02vkbyBMhVQvc7ryR2Tg70WHlirbA99zwLYINfWmmC+PmiobWFQWDDzawCAiF7EzB8nDv0SEf3WsbYsISEh4RrHdaBYbA1jDdx3ENGTRFbEuwDccXzNWh9ETYryTHh6aABUi9lXkYNmzvtjJjKrhlkyTbBZ6COuo+y1XTiOACrTrs1y6XBd9weY6f51kbsMlA5dFE1AG4kyusoFVs+014lJOSm6yGEto/NAJtw+I650NojLiI+pnvNWQZ6pYbRsFyRXXb6E4vb3AwDs7c+hl5U3vlOmkKHRCPSybMVvAG2PLhm7kYm8anG7neF9jCYV37/Epl5j6+IGUixGC4tvAPCbRPROu/1EAM/rP/0UQOQ7H0cDF7tBy3ZkmSDQocsOIaGjtSkkz8xZ1rmsqEtXPpTQblMce3RqnB+qzz1yRRlAKDQaesGmhneDS4fg2tT+sm2sascYN9euRJUOcgnf+Bh1HOsTEJv2iUwpcKZ9MOT8oUuY3vVBAIDdv3MTipt2cfmv7gWAkCLqyJbQBNqV4ftWKhASst1yyeLW/XcIxi6hOOQhdVwwQXk3jrQY6w31aza+4oPsrrcx82LomoSEhISzjvrGkRUrvaGe7tZntcLhjUPnnDa6ZlkyxbLquN0+r6A4riLWKqQBjufhinKujGJ3B5SpYFH7lfew5mxoZdZNeR8d2XhDiqhZCEoGbckZYlc+ICVmu0oE18XXBPUV/RpWrGW0jq9IxTEG69Bcfdd37t+Qdgr+r6CdxlI0q+jHvtQnuqzApV3s6kkfhuIxT4Sa/T4A4NKf/5VfGMm1zV3JtfYr9XmKSmSWVSIXlDvf/e/qK9u435hK3jZuIMVipWbx34noWUBnvnWH/wbgadtq0KYgImQTR/tob19wHRSAV7V9sNmyDAZECdYatTi2zsDilljN92bmukhYjB3c1l2xbRWy6Bqu9aAwGYJ7Hm6dgq7n0+VFFhw/YVtDjG3XP0ZQ9LmEUs+gKAfNuI4h28bgIJll0MKW5b6BLDP0LNca9QPvBQDkj7kL9c2PwfTvmKUYzs8PAkpKBd5RYt0TAGo2bdovUv4PedN5ATRAK3UJVTWQ7vy4wGDo6y3hxxGwSlhcAPA6DAuLB7bXnISEhITrBOstfnTdY5Xr7BNPqB1HB4UzD+cBFcdEyIymXfEEg/SH87Cqm3TOdVkFxnFV5H6hoHw2QTVfDpY5ZvGWdQKNxnrEdBmUnebBSg9qGXpZiUVt8sB7JS43TmvhZqKZnWl2UXSrKKJe+mpwZbz1Vy/cBEOG7L79q95vHxUTz7BVD+UVb8t37bYl1aqKHAf3PWgOvuUPMfvQjwbv3gwAmDzxg3DO9v3DBx4xqyY6T8I6TFMiYysyYdQ2WpQGZ833NFZzbz+b8J7HUnHbwnWWG+pIGOsNdc1jKDeUVJUzSUuN7KBuWVW3Oh4vKyyvHNgyXK4f0ymnN5/3woIy5fldYLVgGDOwAOPdArv299oTomU3FUIfcn+ffmDqplz6OHpXpjvuPGKyou3BZk7of1Z9gmyIhz9J9NU7SCtumiCw57ohF+WYfqJMAcvmfct+sP+u9wLVb2H2NBNmld/5JEx9csG3Y3mlSTHPtQ4SXqsi99+b9OpSRQ7ORO6x5fzI76rTG2oLnmJDYCTNIiEhISFhBOobSFqcHWFhU5S72aWfOau2Z40emSk2s3ES2WwCpbW/rp4ve2kNNcl9que6LANvEd/UNTQBYHjWOTZIyUN6vQx4OQFN4NYq/YuybNCAWyMMeASaWWcpZrdjNb2YQgkb3XP/a9IR6xj7+7AqW+q2Zrt9xuyxXkQuEE45mnVZ+dxQXGscPvAwcmvwzh73gchvf7Q57+L9qOdLT8lK2qmzPbZe54hSL8vO9vTdwxBaQYA9uaJW1bcOUpxFB8isevRPADyJmb+diJ4A4NHM/EfH2roN0fdR6rJCsbfjVWcZpR27x8aQg5LpmOEH6vneFZ1znQ9jzOCybqBSWJ9u7ll4qHQtr+nQtb4B0AyM0vtsVZtIKb9Cm15Wod0j9qQaoJ5irJNOvA+x59gmdYw9PjSwjYH3/OsQGl19r9XPlAIKkcVAhe+oni/BC0unMgO5zaeWTwI7gdIKmIRDiutTeZH7yVfsOht7R20jIr+rD22diuLrI5HgtjD26f0wgI9Gk7XwCoD/fCwtSkhISLgO4DSLMX9nAWNpqI9k5qcT0Z8AADM/QkSTVRedNPqMeX7WD+eJ0ywO7/PT1CqIu5AGOV1WqOaLwHOnPetvz+p0x1rS0iOo1dYBTaIvffqq6+JzfTuiwCyudTNziNIv1PNFuFBOsF6zyR4qjZdZNLusxXOUs2jWOshiOrTgT59W0LV/DHVxFC+ZVUGDm9Tbp30eV94rUqK/qgzQNbKp1bKrRnuksgppWx7wUlMKWY9GqIrc9xuX3kMLKned/n3tgJPNogMlEWUwwhREdAeu4YSLwYcgtlnVQW6beLUwlWWACO6R58UuuLL7ajQ0jCyvjha7X4dm2lQ4xHX0BmfFgVmZdF+sfD4twOS48nyzHdydu6RSJijPCcB8NgkSxOmycbOluj+ILLMeMsDqlOPr2ny2iXVyVq2DIdvTsQ2cscBwP8Wo4Ok/d9wKFgCArsOcYZP2cOLanu9Mw+9J2AC7zvfbI2i008KNZrMY+9R/EMDLATyKiL4DwO8A+M6jVExELyGitxHRm4jo5UR0s93/RCI6JKI32L//cpR6EhISEo4DzEBZ86i/s4CxiQR/koheB+ATABCAz2Lmtx6x7lcBeAEzV0T0YgAvAPCv7bG/ZOanblJorFXIGRMBoKyhgZbzfa8ax9k9ZXoCrjXqZenplVoE2pkAwHC2Wc2bHIvr0EVD6ab7aIoxAVld+11glp8ZZhlq6yTvNAA5U3S/s9mkRaWpSe4z+aoiR74z9XWwbtKZc6Q5Se8ZlWXiubUNrRLHlQZ+FIp+Y3sXxizM1OkUsZYH3BbyH8n6VKM5a2WcQmg684dd3iinYfTN+POdqX+/kvKs50tU80UvbThWcx4MQlwjpuUoSDSUBRHdKjbvB/AyeWyd9VtjMPOvi80/APB5m5bViUDFVv6/UhlImc5O+4fBJSay1NItswloYj8QrVHU2gfiuXMBtHLvV/OlFyoyiaG8pmnicAdfJSBMGdlgHb0cuaCfACs8UATnOLuFWf5U0FDRgEmqGfTznWnLm0baKXzCwSJHPpsGAsn9qjBM9/QNJmNcb7vOXYfaWTd1eszhj8U2Bsyh8gL49yMi3t2hiXH/zi7c1pzvc6YZYSFtgLK+bDaBst+QXs59UKubFIx+vyPut0tgjfVC3BSGhtp6sdcsVmkWr4N5JgTgCQAesb9vBvAuAHdtqR1fDuBnxPZd1ph+GcC3MvNvd11ERM+DXVfj8bdf2FJTEhISEkaAgfoGkharckPdBQDWbvAKZv4Vu/0pAD5xVeFE9GqYtbtjvJCZf9Ge80KYyeRP2mP3AngCMz9ERB8O4BeI6CnMfLmjfXcDuBsAnv7kx3GLgorbYw3YTqXOD+atc7xxdjKD2tkz9WiNrGo8oZTw/umaqWphCJez6VZ7tqR+950Xt7EL3htK5OpxyLXwehJxFHFqkKwoelc9k/VnCL2mpCYHmCy9AHxA2FCAVxfWoRc2pSLWuS5+nkfBUPrtsUbg0W1Xmc8c6jVD923pGuy+BW1ySrkguxaNmxcgtwDY4UErtcyQY0bXfcb32Pq9ou9vm7pknB232DEY+zV+BDN/ldtg5l8lohetuoiZBwUKEX0JgE8D8AlsVz6362Ys7O/XEdFfAvgAAK8d1dIVAgMq8/RScW4P9ULYFzIFldvAoZ29gIaSgkhN8iAYjrUOIqM9osHQ1dGFsap3fG484A4JD7MtqAbd5w0VXis9yDLhGdXXBrlfiVUE5Xld99QkGSygVd0fmbwGDbTJddvC2Ih03ZENYMgOMXZS0dvXOr6Rluu3DbyjSeP5BABY7EPvXwEA1IcHQT1Bny1yIC/89xh6+jVeh33Y5DsZLzS3s8YFAyiTZtHCg0T0rQB+AuYZfRGAh45SMRE9G8ag/feZ+UDsvwPAw8xcE9GTALw/gHcepa6EhISErSPRUJ34QgD/DsZ9FgB+C00096b4IQBTAK8y2UTwB1Z7+TgA305EFYAawFeNNqR7X/CBGaSuwUu70tfOHpSO0lpbioomM9DU5HjC4hDVwWGgTcQ+5Z2L/0SzmyEvniHtoXN7YEYVGpi7fenjY6RrsAqfhQti5IkIYhQUnENv6gvrYSYN573BZ0XuNblstoTePxxFo/UFHK66bmxZR8FYGsql0+89vsb9j6FoAHRr4CKWwp9vvwE+NClyWNfQVy8CMF5NKmvODetraKtVbRvC2DQ2fZpEpxaxRc0i0VAR7GD9ddusmJn/ds/+nwPwc+uXSIGQ6FvzIIDKms6kMlBRgGbGTuHz38AMPvV8GXhKyVX5AABth5BRH0Q/VTDenhGfL+/J/F/9XLhWIG3Th8cHxb15V0gX6Kjr1qAYBGoh73w2Dt51eTIDFfa3pa8GqYRNPX82wLrlDUWid6IYpqxWlTO6D/X0C8qy3gGU8gI0mYFtVDcv5+DFoa9X0p7BRMX1D5tTSnoGxm7qm9zb4P059AmFrbrObq2oax5jEwn+BtBeP5CZ/8HWW5SQkJBwHSBpFt34RvF7BuBzAYzL831SoJHahLxEKcAFG1nDt6NCKMt88JF+5H7UZemppyCWYGDW3FnnBjPiUUa6rrgSX24/FRWU5ygDpZCpEiwphB46i4TnjIe9jnXd+3w8VeSedzEJy81UN71xTEbqbVJQG5UVPSeZo2y9uvtm03F/Ue3vxRmjpfcTDDXrPKC4XPrrXICm93CTBu2iAJelzzPVq4Gsex9D9xTdX2/Z2whiBABm6GSzCMHMr4t2/S4RveYY2nME0NpcJGvdeH3kBShvvlheHEIfGK+P5aUrIKV8ZLJMpOdrX8ttc/sfQucHsEItb6cBd4NFbYJrquY66gi8WlnHIgx6hAzk07Up13LipBS048VthPiYwSV+lr0Cbh1sidM+CoYWvR+NsQOmykI7RdfEQP63v1VehHYulYkA2AxA2QrQBADKR+QgHfk9DQ78x/wekzdUB6JIbgXgw9EdP5GQkJBwQ4CR0n10QUZyVwD+CsBXHFejNgJh3Ewi9n6SQWRVCa7s2tqH+ygvXza/a418Z4p811ImYw1pMVap3keZJfXSSyO8YPy5zoPFXiOD/yYiN9BAbqSgPudN5lCFKwzSdBaUC200ORm4B4ycicp70NFstw9HeR/XCwbvX9KJKtAeAg++yazRCmOvPaFpQCnvpBAcd/9zF7y3WV8+8rnbBiPRUB34O8wchDsT0fQY2rMxiKihkXT3kqesdatz+YFPa+Pl4ZeWnIfeT7NJOGgN2QViHJV2wgiaayxXPVBmkO8pb2ipzutGDEIUvwelALdWgrMR2Xfm3JkB60a7zuAi76Gum3ej9Xoeclt4T9cqRvWfLpuFs53Yd+GowuaaSEAEgloFv4Pzjtrfu9p/wjCaxfHXQ0QvAfDpAJYA/hLAlzHzRXvsBTAT9xrAv2DmVx5XO8a+kd/r2Pf722xIQkJCwvWGE1op71UAPoSZPwzAn8Nk6AYRfTCALwDwFADPBvDDdt2hY8GqrLOPBvBYADtE9DQ0drebAOweV6M2AwVqdJd24Rof0CjxeWJb5juifNI9S12DktpYOzhq2SNnXZKGYt3h5TR4bUf9LgeQpzQyr52RUkEsS1BWPjEG9T76qaWxCDolB+A8d2SZA7EEXfvXyaN0XSBu58AzDIzYezPzzrq+p+kscAqBygINkYCGtoo1l6G2rWr7AIa05W2DmbE8gUW4BzJ0fyaAn7Ypkv6KiN4B4Jk4pon8KhrqHwL4UgCPA/B9Yv8VAN9yHA3aGI6GEu6fQHdnGRQa3tNDXJAX/e6nfR9A17lD+1zbxgxSRxEqIygdruthQdFD83XWF7nZStdKaB3aMJzHlcrMu+yIAjbHVwhGR1+JaH05WDXtWWPQOsKAdWoYbHMxLHRdksLJzLjBLtp2ICfwg/v1btNRXq5IqAy1bxsTn14qdItgrJXu43Yikvnt7raJUNfFl6PJ0P1YGOHh8B6771iwKuvsSwG8lIg+10ZWJyQkJCTArJS3hrB4kJmf0XdwwwzdXXO6Y7OirKKhvoiZfwLAE4noX7Zaxfx9HZedDohAxSSYqQLRzEKHqbUJEVUBCG+PaCYk/dAl1tAcNp4xHUVT2IBOWT0zjmIuRqr6XNft9orU7+ijqNaZ1ecFyKaW11cuBgZ1X3ZXuau0jhjr0k8nTFeN1SA9bNyLh4uJUJl5F7VbFS8LtITWc3LvsFoa90kZdzHmGa/zPY3q+6KvHgMlta1Egptk6IbRJB4vTnscgPdtpUEdWEVD2URJONdx7NryGSMydFGrAxZhJ6nrYMF5f0zX/Wp5LCiO8uEfqeNvXn6rvi1RLhwPMl1wy2923YMTClobe4OtP6AtesrrghMUAIwgioXEGApxTerQ172Oq+1JUFXreBzpSFC4a/MCVDSOjzTbBcllVSO7Ejl7kcpC2nHo+4mFSN9zXNfWIcsU5/KWqCkGn0jW2b4M3QBeAeCniOj7ADwGJkP3Hx1XO1bRUD9qf76amX9XHiOijz2uRiUkJCRc8zi5FOWdGbqZ+S1E9LMA/gyGnvpnzHxsRpqxcRb/D4Cnj9h3inBxFkMzUm2oJ6FZkNA6Apokzo8zchbTOWsZ8jxZAxulrzjCDHbMbGyQrnHPtqMcrzk4A7Sum2BAlQU++fK9mDiK7pkozYyDnlucB/kkpCH7aJBVBtfOVCpbcESI6zlOqqqrvTp6L8LJw7VF7eyBswmILC01OfTvhuvaOCWIwDte9nuc+Xe6IR3bDoZVg+le+q6ndRO69UAzsKxOxBuqM0O3PfYdAL7j2BuB1TaLjwbwMQDuiGwWNwG4tnwGlQqjgQcQWIX6eMyuRGsjENtMBuvoqLPzmm15dBxhMFprIIs8yxzkynyAjcy298ll2Qw6WRYIa4p/R1C7N5lj0xnqS82aXCbyeMSzU1lHUr01vNtcfSu91LacP+yokHVo1aaiYL2assxzzv7doHEp94Jf0Ic6Lxp7EWDoLOeldtR767KtuPfnvr+ub+4YqL+0+FGDCYy9IgdwXuy/jMbXNyEhIeGGw0nZLK4VrLJZvAbAa4jo/2Xmvz6hNm0Est5Q2/I62dRzaZt+3dxjfO/UXgawzRxHo+qOn51Pt52FM31dN7RUEXqf9a7wpzKroViaZO+81yj1/hVDLU5lvqk46FL7+whjZ/oNsN2rra2epa7WNE6WVuyFpWf7ghM1KRDbxY8EDUjKrmcfLKJkVzuczsy5noocqamPeSYdmmuT6ucYNPIeMANVEhYtHNj8JE+BWc8CwDW2+BFZGoqOyctElMsDdZCKVo3bsD3EuqHLWJs/LVTssR/CJgPSCN53rXICe4PI26Rr8ZFnwUA4WI9WUM7rae8W8MKmNl/OTRBZn61F60Y4rbJJjUyJHl4zdGwNCmpT4b7pREnX4FqFKeSd8M6NJxSxFbJaB3SStFmEXoPFZmnW1xEW8nxSZgyQ57nv5hhxI2kWY3vwTwJ4G4C7APx7APcA+ONjalNCQkLCNQ8XlDfm7yxgrGZxGzP/NyL6OkFNXVuLHxGB81mw3cJRE3p1lRlX0WX3H+MxNFQOM8Daay0kZ0xDGsbGlNxw9t4Yvfl3WnROHWgbct3v8Lrh5xX479fLRluw3lDSc4qioKxBmnB0KokjrMK2EfV0zMZulYXvBvBr0evpHkCq0ZBdnjT7O4iJUVmjgW9Jw+/SzGPt3Z9Hqv2Num/+GDQMxsnkhrpWMFZYOLeGe4noU2GiBB93PE3aDAwCy0Curs66Zochcf4gndR3zF0/dLxvmVBS/npiDWgCs/LX+bapsa9wBX3W9Wxc2UPPjfVqakrk65Jpz+MBql1/R7m6Bte1X/JWRh43rpl2Ow+v62zTulhn4B4dIHnEgXUbwkQrsFgomTPrHpvPgndM0x2Tyh/meQc5vNyADTF4b4qBiVnnhAww/UvWGX9DW8aa6T6ue4x9m/+BiC4A+Fcw63H/VwBff5SKiejbiOi9RPQG+/cccewFRPQOIno7Ef3Do9STkJCQcFxINFQEZv5l+/MSgI8HACL6+i3U//3M/D1yR5Sj/TEAXk1EH7AyMpHIzISj2Yycaaw1u2A9Lp/JJrOnVRpHVDazBkhoE5r6Z1cbokV79Z7YkcV31XOV70BSCDrK3eOrGCjPUSbB9piYiBGBWNv2ntl0xn9cThpDUKEmxjJQUbwOUsqvgEhKAVmozbOkoYg28t5qacDrPI9Ys3BlAlunokzW2URDjcG/BPADW2qHxIY52smozrH6GncYYFyn2ebaunF9Lhp2oB1MqjkuhQYAxJztKhpsXQyq7lnLzuLpsZHl+rJV3i18MvS3XdcYcnNtPYt1BpqxlF5f205jkI+wjvdd8H4lbST3xx5F0q3ZXtNbp4qEx6hGtc/vKz/un51C5jhsl6Kcs6I1jMFRevdGnnERvpaI3kREP05Et9h9jwXwbnFOb452InoeEb2WiF774MMPb6E5CQkJCeOgGVhUetTfWcBRNIuVInUoRzuAHwHwIlvOiwB8L8zCHqNztNvFQ+4GgA9/6ocxZ0W32jt2ViNnKvHiLWtSWOHF0ezXzmxWzQJHzxLl7Klr1tTT9sF7kjRcfJ5q7xv3fBqthJF5L68+UDzDBfodAoBBCvLYZv2srwmNYhMEz0fG9QCgamF+OGrO9asuI3YXrdql4Q88p9a7ctd3aBg+5qNVSKQpyXKV6vU83BRrLn503WNVbqgr6B6oCcDOqsJX5WgX9fwYAGcX2SxHOymgMK6zKz2XxgyeKhyoVnKeut1Re8vv08nW4P6De+yJmO6t35XRV+eAu2FLKLgPd6wwdU2Nqa6ONgf01gA11stxywHriB5H/bavYbvExjazsdjWAEhZ4G3nd1dzIJs0noYt6orakyGgoaC6+mzrfYn3RAqcCdtjh5DpncTEkwt5LantZ7O7wbyhVqX7OD90/CggojuZ+V67+dkA3mx/vwInmKM9ISEhYROk3FAnh+8moqfCTBTuAfCVALBxjnYicFaMpASk54+IpYgNtfFstnMWPkCPtGbyKzSAnlloa+asohlTjD5791D7O9rVdd8rn9E6YC2M9ivK1NTU30E3NNuNB85KGm9do6tr3rirWgiuOwqtKTFEyw1g8H1pCuhMJgUUu811cb9VHdoD0J7Z+/3UNowrFVzDztmgh5IylUlNYkXwXZ+x+4hIwuIEwMxfPHBsgxzt5AOJxjci5Glb9BCLAcqe70/tc38dGsjijhx/6HJQ7hrxO1TzIwU/dd0PYCg4OZAD/QLEUhCDHlFDrsJ9bZD1s24u1bpNe8SD1RgPnBHPrFPYHNU+sa5wPapNpKee1r3J87LGE49YA/WyCdLLJgBV/tigXSjup9aLilXeP+FxgsIO7Czc4XlgsKfI/tXpKbVl29KNFpR3mppFQkJCwnULZqA+I55OY3BmhAUDqBhQK1TNwJaGAR9xMaPzuZjEjLaPjumjkuJjwTV9x6NjrYAncc1Yz6mWcdZ7lmSBKh/kn2LttSA/29fRLG7IaN9lAPXtzfx5PjDQ5cIS9ft6qeOZxTPTrhnrkHbTalvHzRzRuQBw/SjSVle1bYQGMqilDN1vnxYsNQbbD6hetq7hOAg2fkain3ptAvDBs4P92dJQTOS7pdY8SAESlP/+WyvoHYsCwODjKfiaxJkSFjUbo5NE/NkHSYyHBAtlRrWF7firXEW57b46qOoDrQ+5N8jIe510CAZSwQfVuo0WtRYNUoE7oqtfmacYC0u3TwgPoMMWMsIzzFTTPdC6QTUWHuGNdXg8DQqCaPBYh5JQ0YA4RGn0BRQ6j7q6oXBEY/rr7hK00cteOw1+j3AaTa3GFGjPO2BSQGaHGJULO0Rks7B9WN6adt8QMxzT4373uWcqAtxRIrLbrryOi7YATjRUQkJCQsIg2Gg7NwrOlLBgZtRsZhUOKpq0aTkrG3jPighstRAitGd4fTQUurSO0Bjchxa1JVRzVnlr9gWY+5GzL9N2UUSwP7wnEkZ0htTCLHXUo1m1guVahs2RPG7Ly8sW6Wm/nhlzPLOP27NOnSPP7Q0a60KPZw6Je2KprXVhFZ21CVb0w9FpW5yhmtSwoToyart8U7H2UOtQO9U9n6hmq130NUv8Jmp/+9v2hWKMYgnPDM6MsGA0nUyJLqY7ukjciSCvtcfqSKi06KxosAg40g1yJXVCeIFIzwv30ch2B8UGTBsF9+QOKRjB4eWRtN84e4E7eche03EPYwezlodKHJw15NUkvXW6XChXYV03yiHBEbWruf8saFvs9ts7qQBawYvBaZuMUHq1V1Xfe2vZaIRba+h9Ro2tQWWmz4p+ytYY7KikgF6K6+x4jfJ832bxHqV/vSIOYzLFuV3f/6ZINouEhISEhGEwJ2+o6xWOStHczB66JH89MBnoO0bRTFQa0oVdLaB6YuioLWoFLeWv04Zec5pFPCuL2+n2mHaxn3Epau6jZjbbot2+DJWB0Rj4TUV9M+KsNfttpU/vmdHKWSkTte/PszncohYdhaZU1k5DtOXZXuzfH9OATbtM6gu/UFDcjkDriDSEnmdtyukZkNbRMEZQg71PTWpV0gMqMlRXUKjsOywr3dIYGvrJbVutY6hNq15loGWEu6VWrQlQvr7tqBaGhkqaxXULKTA2udah1Z3s4OoghUqgCsvI18iWEEPRuEbGanvNDLeaYywsMkGndw3RgWpOzX0atZ1a+812mP9n1cp4fYMp0Oag69oJCA7uL4akFuUgQAjfSyzUt0I5iPbE75M7BqumbbIdoSABhDBpBeqFGQaOstZIF8bQhEOJ/WIvJifkFzXjsDRllyvcXCXW/Va73ikh9HT0wuEYorY9uN2/zzLOnLBISEhIOCkkzeJ6BLeN1A7bep9SmwjqELOL2lJG8rc08ikCMhc4FLUz9OZoz4jcLDYuU7aHNAeUkubuGZ7xSydk9tQ8Iz97jw2Bq2bvUoOqOWynhHwPNTNYzMw4Oh7Uh4YmI8v5uWcV+9OT0NakUbOr3FXw9v2Bc8I2c3Bd2C67zWHZ5jxLW3VobLH2MdjeoXOCVDIr6M+INuyjC6uaUTNjXpntg1JjblVeFenmQxP8Ve8kvpa5/f0oMn0pdjxpnF58ejFPR20DSVhch2A0A2if3WFVH1m1LERfv3Auu4BRvz1vqxnLiv2AeFAa60GRucGYvODw2/Y7zVRzzO2TnK+7V7/MQOD50XxMoUdKbDNpvKEKRR3Cq6k/G/jaa2aUfgBp6ltUGqXWKDrSg5daoxQvShH551IoZeq0o4gCeSFg7C7yHnlQsMT0VVNf7+1szFzUYmCK64wFSdiGtk0maE+rVHfegEdehHXsONJG4wR6Q3tqzO2NXpzXuLKoMa8c9aR9PykyCt677Nuh4G9/r0OCRhGgXV/gps+6iZjr4wrUtlM4ubslZoqZUdfJwJ2QkJCQsAIpzuI6BDOwqLhzVjjovSIQaxbyXEmTMANiHTnUWmg1moMZ88OHJR6ZlwCMZpER4VF7k876tdBQMgKmuZmKnZ/kwcxMXntQ1ihr9rPyq8sat+8W/pgiwoMHZeveXVvd7P3CNMduYWamV5YVLs0rv31+muH81HSVQhHuu7r09bm2uKUjrywbb/ey1ig1+3IKRf53qTUuzStcXpj0F9NcoVBuVmp+F/aGze9GG3Oah9t2jzGzlJTXQkQciTHahzRFnwaxZvKMFkItwVF7/XWY82W/FMfitgU7QtqwfXxFO1sluQ1BK1oqVWrLlxbmHb/n0hylbtZ0yBQ171ArFKp53pkClOWBMmrepYPs22GckNEgHDSaZ6eJGy9EUOh5BQaxjL1oU1RHBXOK4L4uUTPj6lK3IjeDaNBokI8HE4mY6inrhk7q8tZxFJPm5netGffvL/11h2WNWjMuWuFx687ECwTHt1ZWrZ3XDX1z53lgt1CYWq45U/AC6cGD0g/Krpz7ry58GbIzt20IjKWtTwoO89ErP0DPDhRy+zVnZOglV1aXN4invbShp+4X50hqohRRyvOq2Z8RoAQNVyjy9Tsh0giWhu5wv127iVgMTI6Wgj3WtHebkb7KcX+2ji5qytfjG8HtNgwIl3YZ3YKmdW7/IdsKA9n3jV0C2LdeTu++NPeTn0WlUdbaC/VA4FvbWSH6lGubEeohLWWL9wLfCw9uhIe71S6bBQvBIRF+/9sf2JPNIiEhISFhGJyExXWJUjPuvbrw3kaxR1B8LhB6JsU0DwBv1HPG2KvLyh+b5maWv6hqlJo9DXNQ1ji0msXhssbhssbSHqs0I1eEidUmHri88LN5Ux/jytzUcVjWeNIdewCA23cL1JlCMTHn3nPxEJfseYuqxn5p6nF1uvoWlcay0r49y6r2lEGtGZNcYWLvY5IrTJz2oCho584kC7SOWrOvY1lp1OL+a91eatJduzvJfJmm7uaBH6JuPQuHm2YFNFf+PRVZt9YxyxWmmcIslzNdd092dmvLbMWbdDgIjIW81gU7AoBME9MERra9KKQGYurvqqOvTbHTwuo2dpbCzTOX3na1Bg4r7fu+5uZdzysNzYy5zaQ7r5pnXyiNPFPe267ImvThGdntQJswx0pwB2XlvlfzjdZiuw/OKD/EBhwd3KlZn1WcHWFRa9x7ZYFpbjqo48ZdhwbMvpqBhw5Mbv6dovEeych0PtfZFZH/QA7K0IqlmfHIoVHFS804FIP11XmFi/bYpYMlrswrP7ACZsByg+TuJMPOxLyCSa6CQThThEvW1vDey3Ocm+a456K5j0uLEg9fNfdw/+UFHt5f4KI993BeQTtqqTLeGtoJq7JGZe+lsu3NJ1ZY7OSYWVvHZJojz1UgICSkQKgqjbrWPu2BrpsUCM5TJHMU0jTDrrXX3HZuglv3ptidCGElhIcUFlfnVcCL33HTNGiPFBy7Rebf6yxXLbuHEx4xDRIKC4MxFHfgDht590ih0xxrJioekTtt9xrCqwcleYYrq89Go4BWdLUUEA/Yb8TZoxy16uxggPlmNKPp+7rCzsQ9+wyzvLETOFoKMO+iFNvu2wOsQNcDdo26W0jIpITBtv1f1g3NWm7Jg4kZ/tu6EXBmhEVCQkLCSSMZuK9DMAxdtJjXqBkotTHyylmwM+K6GcbFeempFyAMUpOoGZhXIdXj6KJLhyUuHixx1W4vlzVK6y1SlTXqWge8JinyM+28yJBPhGE3U5jamdnNuwV2Jrv+2COHJS5ZjeX+ywu895ED8/viHIdXl1jYY1XZaBJuxu80jWq5gC7NjLFeHkJXS9GuDMXeBQDAzvnzmO1NoDoeBjO33AW15kh7Mfdfzg+gy6X3/1fFBLM9c0/3nZ9g76aZ12bO7xY4P3NaltEOpKYhcSg8ruT7vPXcxL57G+dRG1oKkB5WVtshGhXv4rYlQtqpW4MArNHV7Y/uYSg/GdCmwnSPZhF6+nWU0+GN1UWzaTSxFAA8zTnNFcq6iaOJ221Sz1gHjqWgOScMzZnQLBqtrrIODM74bZ63/a3Db3AdDyYZlOeM8Fett12syWwLJ5l1loi+EcBLANzBzA/afS8A8BUwCum/YOZXHlf9Z0ZYlJXGex459FTOgR1QlrX21I6kMtx//1t02K7rllXtudrDZY2FLX9xWGExLzHft7TUomplonQdygkN5WwDucJszwyWO+em2D038W04Nyv8YPnIYYmr8woPWerp3kuHePCRQwDA/uUF5vulH6ClcDCD+BKsradW2fw2x2svMHS5xPzSAwCAqypDNt1BMTtn2jnZAWVhbqi+ADDWta+Da1uvraOaX8XhI38DALikMkzO34qd8+cBALO9CYqpdbGd5phMMxTWXXdnknlBsjPJg/cGwAv8C7uFCQwTPI6b+ZWaUShCVjdUiKpCT50++1WXMHHIVOTqKZ9FR/SXHPSbIMvxA44iagLPLH22tBHU91ovOIeZ8E6K29x1T847TdJOpWaUtfZ9v9QcBOE5uxzQ2K8cnF3MnEvIbDJD59FWShpKCo4VwiIT5y4q2+/ZehtGAbByohiPA0eFmTidjLAgoscD+CQA7xL7PhjAFwB4CoDHAHg1EX0AM3czmUfEUV3KNwYR/QwRvcH+3UNEb7D7n0hEh+LYfzmtNiYkJCQMQWse9bcFfD+Ab0ZomvpMAD/NzAtm/isA7wDwzG1U1oVT0yyY+fPdbyL6XgCXxOG/ZOanrlPewbLGG9990c/6a2/k1Z42YQ5fnFIUeJnIBeFY5KqvytqU4zWEZsaqrQFdB+daz5FqCV0u/cya7Ix9smu8nIppjunMaBZ5oZDnCredM8bbSaZw/2UzU1xWGlcXFS5ao+PFKwuhydSBgbkqQ21BahZy1h+DdY16Ofe/y4PLWOaPAACyyQzZdMf8zicrs87K47LOQJMBoC89gOrwKgDgYLLj68gnBaazwmtdi5umWFbGML4zqbErvLMmeYZHnZ/65zTJFUoXlKdYWHE1NJP3siljo6qY0cr9ZrttcPXGcR3O0OU5XUZpqUW4GfqlRdXKa1QoqQmEMQpu1rxbZJjmKvDEk/U8csheW3AxL4DRmuXsOlOESaZ8zE+hGi+mg9L071LE/yyFhn3poPT1Oy86BzmbnwiHiTpXbQO30PBb2kWkdWRWf7u6rHCf1aacdimDY5eV9vTwlXnlKalD4dV4VPR9T9sEEX0GgPcy8xsjr7jHAvgDsf0eu+9YcOo0FJm7fy6Af3CUcqpa48FHDlFXjKqsvXooB1LHt3dJetZGpdSCMmo8fDTqqgoH3brpJJKi4VoMjpUZrLXdzvKJHxB9HY4yqhnzwxLvtnXmkWdQVdaollb9X1RY2o5fV8ZG4T2gloeoF4fNbyuwfLujzt20tfS/q8VhcB6pDKqY+HswNJTyx+T/uGzz1/YYySYzqHwClZtyJ3vnMdlpvLEmOzkmlobKi0Y43HZuisfdsuPpDRlMGNs4tGYR7gsoBjIrSLRwY2rcNM2pbWERChLnsWOuoUCAzAUFeVVw+JcXlcmhZUn/eRXaXWLbWh+ki7NzI5aBjtp7/5hATRd0uay0p4vcQOo92uz/qXCVds/SD/DCLuHKubKovIAGgAlUr9t0XFZMJWY9wqLIFAp3kjJZGv7m4QN/jbvOBbxKQSYFxNV56e2Mi+V2BnhmHdj9VuB2Inqt2L6bme92G0T0agCP7rjuhQC+BcAndxzrsr4cGy926sICwN8DcB8z/4XYdxcR/QmAywC+lZl/u+tCInoegOcBwPSW9zv2hiYkJCR4MIJJ4wo8yMzP6C2K+RO79hPRhwK4C4DTKh4H4PVE9EwYTeLx4vTHAXjf2Aati2MVFkPSkpl/0f7+QgAvE8fuBfAEZn6IiD4cwC8Q0VOY+XJciJXMdwPA3mM+gA+uLjHfL7E8bAy+zsgLNLPortmwnwWLly81CflfagvxdQHtUhmqSOVmbqSnO2aWbmfTpbCiOg8q/+w6rEm1NWQaLUN4XC0PPZ1TCS+n2moIrh1DKrOONAl3n+5/vTTaShmd4+5HalcOTsvqqjeb7oB1jXxiNC2VKUx3THec7U4w2yuQ23iJc7Pcz3rPTUMDd2cMiHKWYwSzUttyc17NfvbqjM5aRNPFhmQtKCln9HXw3lc2y67fL1LEOOonNgDL//G9xOcBaMXsSCwrHQRLLoKATI1adwdOuuDN2GvQ1RFrbLIOFpz8ZJLhzpvN+3zU+WnwnuI6u7QKoIm3cR6DZa3h2LX7riyC8+M2S8cWE/tU+e2r8wpL+zv+1jYHHzsNxcx/CuBRbpuI7gHwDGZ+kIheAeCniOj7YAzc7w/gj46rLccqLPqkpQMR5QA+B8CHi2sWABb29+uI6C8BfACA13YWYlHXGlcePrQDadkIiLoOBj3WdTAw9lFLqmi4+cwO7p3ce7lEvZyjsoNpvQhdUlU+QbFzzm+XwpNIVztYLuxgGX34WnPQNunJxKL+uloGdeqquXcd3VOMrgHe7aesaWc84LOuA0GrqyXQQQMP2TbccyvnRsgtDw897VRMzSDUJRDc4HDvxbk/5qiGnSKzgY7tQL9JrlqR/UF7ogSTGTVC4G8uzX0552Z5MPDtTLJWCnYZEBZ61IWDW9dvoKGF5P26Z5D30Dfu3IUYyA/LGsuqERZd5Tu4oFB3rgy6XCxrLxCkq6gTFK7vLhaMdz+4DwB4z8MHODfLvQ3u3DT3AZgxYqEi27usdEAzymclhWOcreDiwRIPX15gfuBcypsJ1taEBR+/sBiunt9CRD8L4M9gvsB/dlyeUMDp01CfCOBtzPwet4OI7gDwMDPXRPQkGGn5ztNqYEJCQkIfTlpYMPMTo+3vAPAdJ1H3aQuLL0BIQQHAxwH4diKqYAJNvoqZHx5bYF5kyIsMrM2MRmsO4hykoVp6Kg0Zf502Es/0AWHELsMAN/ffUU7uunpxCHdmtjxsAta84TikxlwdVWS4Luf7/h66Zv4SwbrPom1SQ4i1DKUyQGoPG3wUKtIs4nZkk1mjZZVLb7TPD00wlxazW2fwvzKv8M4HrrbSjQBAtdTBYjRZppDllj7KVKC9xf2i5SWnyNNgcZyH02CAUHuJjbZAm2rq0ywWUuvgdn4tf08q9L7qy6flvIJc7q9MSYeJ7ngDtz3JVTCzrwJPwOZZOc9BR+NlqnlnCob6cd5I7lqgcTBx74MUeW3w/CzHuVnhU7ZMo7Qz8rkdLJs8bM6AfWjrKxc1ykXVxB/V2tO424qNYOZ1DNzXPU5VWDDzl3bs+zkAP7duWQTT6ZxXU+c5ipBPCrC2tz2Zio4f8uuxVxOqpffCjD2FnGcPeo4H92cFBmAGSDdI99lRun5nk52IItIB9SQF19iZT+zdNVY4KNXfflIqEJiUNb/zyQ5UPkFuKbpiOsHM5o0qphlUrlrUHNDB60s+fQJQ1QxiZjVQKyysAOgbXOM+o4U33IFm7wGUqWXgBioHs0meBccmmQrqk26lVUxDcTvTQNc915p9Ir22ANKtfY5yiumrLgzZgWTCKhnBz8xwxIeOKM/YLd27wFrBnXcI2UWlgXkZ5EiL2+KotYXw8HLu8s5j8GRWsOONJlHXK05bs0hISEi4PnHKNouTxpkRFswm/iBWdYdeZhft1OXV1EXruLiDVpkD3lRd9Q4FyY1tNynlZ+wKAHcY5MfA3/sKd8CYspL0mdMmAHgvKe/xlE+Qud/FBPlk6lN6THcKTKw31HRWoJhmOGfzRsnsvE2K9GZ2GVAmWUNFyVgZN+OUXmaSBulPAW77kutPilBV2mssh5EGIL2H5Gx+kme92kKMIS8it8/9l2nnKx0a1IfiNeK2dFFo8jhb7ySz/nUT6EeRq3/8TN12HmhgKgjg69IertjcTrFBu27FPzXxVFpQkibtTXOujGmK0/FsCsbJ2yxOE2dGWOiqxOEjDwIYNxiPQZxHqfe8nsF1KGCtq22rt1d38i4KLOuhtjZRoSnLoALh0AToue2YdoqFh/8d2xA6PG5M+ykYaADjXeMgo3bb7qnh4Fnp8Fz3u+u/g6SovG3DnYtmgKxgBklJoTTUTt0aINcRJI0QaHsALTsGRCMoOXiWsa3AJYpUiqAyhazH9pIpQmZtCua5hd5fnbnW7O+u+zXboQ0ltt8EEwBrd6iWOkhUWVfap9zX9jyfWcGuEimpZDkR3AqY14mzuO5xZoRFQkJCwkkjaRbXIXRdYX7pgV7qpw9dnkIAWllWJdahiLqu3wRSq1ilgTjEmoM02q8DUllv1llHOwWahjDaq0GDfThD9Sla7CzxYN5QMS64auqCtlQzC/axFJlCVoSz8iG6pytADgg1la7jXdt99fVpDHE+KWnglnEOLtuxTGMR5AGLFptyOcqct1dg7CehSeTNynUqV8gyblLZqyaNPqmQKpLPO9YW8hXPOowlafpCHH9SCW3Cpe9x9yvT3ricbf5ZRKlt4rxodRSce2Ssl+7juseZERYg8oNUHw3VTeOU/hdFwVXxID806I8RIMPntNu2jlAKAg0jdXuM7QRoC86u+1VCSMjznA1HXpdF7sCxp5dxMjPPf3GYeRokLzKfphwwCRedG2uWG5fWXATbxTaCZmDrpkvc73jw9rSQUsHAN9bW4NAlZJz3kxcCtaRa6hb14qKNHS/vc4hFvLwUCEPZTRXIJMcCoDWASoPFvbBmMNt3SoTaJr9yXmS1EDQOFD2LLrquz+VWJvWUHowuKWdDJzXnxUk7pfuzaU/TD3W5hMonAQ3lgmvrLQ3wxmaRVspLSEhISBhC8oa6PmHiLJoAsy6vpL7tMRhK7x2eNzzTWJfCArrppL7rhoz7q9o/6rlEGoebzXlvLGngHtA6DEVFfnYqPWcARyvYmfaiDoyyjkYBTGoO/9t64DiPJ2e49cc6DLdd/+PjXfTKmDxOcdqOeAGe2KtHag8uZYjTJJz2UFc6mIUD/UFmcuavSDxrRfZ5NGUqRSgXzXWx1tAF1mH9UrOJMzzH8U86Wjgo1kritCKy/CBFT0dclDwWB9aa+zs6LSzruVFwZoSFc2Mb450gA9p6y1uTMhp/7TiBEFwzcE9jorX70KLZuqK4O8qU7rGufdrGpbNqAvqc15RLQthFSXVFkMcUmPytooFMurxKt035H2jSjnceEwOnKzMOCFxn8OzOoxR5fOkwHX5A03TkX2q2+/tC/D77KCPqoJSGEAuBMW3roj7HTnK6tvvKbNrYf/5QmUcCp6C8hISEhIQVYOYgW8JZx9kRFvbFrTt72Mw7SA1qF0NxFX3HsiEt4Bp6S41qX64812Gbav9REDswhMf623ga7R9sT0+24FXXnQY2onw3iF04HToo2SyuS5igvPtOuxkJCQk3EJKwSEhISEgYxg3mDUVxaoXrFUT0AIB9AA+edlsi3I7UprG4FtuV2jQO11ub/hYz33GUwono12wdY/AgMz/7KPWdNs6MsAAAInrt0Dq3p4HUpvG4FtuV2jQOqU1nH/0Wv4SEhISEBIskLBISEhISVuKsCYu7T7sBHUhtGo9rsV2pTeOQ2nTGcaZsFgkJCQkJx4OzplkkJCQkJBwDkrBISEhISFiJMyEsiOjZRPR2InoHET3/lNtyDxH9KRG9gYhea/fdSkSvIqK/sP9vOeY2/DgR3U9Ebxb7ettARC+wz+7tRPQPT7BN30ZE77XP6g1E9JwTbtPjieg3iOitRPQWIvo6u//UntVAm07tWRHRjIj+iIjeaNv07+3+0+5Tfe061X51ZsHM1/UfgAzAXwJ4EoAJgDcC+OBTbM89AG6P9n03gOfb388H8OJjbsPHAXg6gDevagOAD7bPbArgLvsssxNq07cB+MaOc0+qTXcCeLr9fR7An9u6T+1ZDbTp1J4VzAoA5+zvAsAfAvioa6BP9bXrVPvVWf07C5rFMwG8g5nfycxLAD8N4DNPuU0xPhPAS+3vlwL4rOOsjJl/C8DDI9vwmQB+mpkXzPxXAN4B80xPok19OKk23cvMr7e/rwB4K4DH4hSf1UCb+nASbWJmvmo3C/vHOP0+1deuPpxIu84qzoKweCyAd4vt92D44zpuMIBfJ6LXEdHz7L73Y+Z7ATMYAHjUKbSrrw2n/fy+lojeZGkqR2OceJuI6IkAngYzO70mnlXUJuAUnxURZUT0BgD3A3gVM18Tz6mnXcA10q/OEs6CsOhaweU0/YE/lpmfDuBTAPwzIvq4U2zLGJzm8/sRAE8G8FQA9wL43tNoExGdA/BzAL6emS8Pndqx71ja1dGmU31WzFwz81MBPA7AM4noQwZOP7Hn1NOua6JfnTWcBWHxHgCPF9uPA/C+U2oLmPl99v/9AF4Oo+beR0R3AoD9f/8pNK2vDaf2/Jj5PvuxawA/hoYSOLE2EVEBMyj/JDP/vN19qs+qq03XwrOy7bgI4DcBPBvXUJ+S7bpWntVZw1kQFn8M4P2J6C4imgD4AgCvOI2GENEeEZ13vwF8MoA32/Z8iT3tSwD84ik0r68NrwDwBUQ0JaK7ALw/gD86iQa5gcbis2Ge1Ym1iYgIwH8D8FZm/j5x6NSeVV+bTvNZEdEdRHSz/b0D4BMBvA2n3Kf62nXa/erM4rQt7Nv4A/AcGK+RvwTwwlNsx5NgvC3eCOAtri0AbgPwfwD8hf1/6zG342Uw6ncJM5v6iqE2AHihfXZvB/ApJ9im/wHgTwG8CeZDvvOE2/R/wdAQbwLwBvv3nNN8VgNtOrVnBeDDAPyJrfvNAP7tqn59Qu+vr12n2q/O6l9K95GQkJCQsBJngYZKSEhISDhmJGGRkJCQkLASSVgkJCQkJKxEEhYJCQkJCSuRhEVCQkJCwkokYZGQkJCQsBJJWNwAIKKrq886Uvm/QkQ327+v2eD6ZxHRL695/iUi+pWe4/8vEX3euu24HmGfxceI7W8goncR0Q+dZrsSzh6SsEg4Mpj5OWzSLdwMYG1hsSF+m5mfs/q0zUFE+XGWvyU8C4AXFsz8/QD+7am1JuHMIgmLGxRE9FQi+gObmfPlLjMnEf0mEb3YLirz50T09+z+XSL6WXv+zxDRHxLRM+yxe4jodgDfBeDJdsGZl8QaAxH9EBF9qf39bCJ6GxH9DoDPEefs2Uyhf0xEf0JEK9PNk8EPEdGfEdH/hsjqS0QfTkSvsVmAXylyGX2EvZfft219s93/pUT0P4nol2CyB3e2h0y205fY/W8ioq+0++8kot+yz+DN7vn1tPuTbf2vt3Wes/v/rS33zUR0t00BAiL6F/Ye30REP00mK+1XAfgGW19vXQkJR8Zph5Cnv+P/A3C1Y9+bAPx9+/vbAfyA/f2bAL7X/n4OgFfb398I4Eft7w8BUAF4ht2+B8DtAJ6IcHGjZwH4ZbH9QwC+FMAMJlX0+8NkAv1Zdx6A7wTwRfb3zTBpXPaitsflfg6AV8EshPUYABcBfB7M+ga/B+AOe97nA/hx+/vNAD7G/v4u127bvvfApq7oaw+A5wH4Vrt/CuC1MAvq/Cs0aV4yAOd73sntAH7L3RuAf40mXYVMm/E/AHy6/f0+AFPXFvv/2xAt9GPv4YdOu9+lv7P1dz2o2QlbBhFdgBlsXmN3vRTA/xSnuOyrr4MRAIDJWfSfAICZ30xEbzpCEz4IwF8x81/Y9vwEzOALmOSLn0FE32i3ZwCeALMIUB8+DsDLmLkG8D4i+v/Z/R8II9heZSfnGYB7bfK588z8e/a8nwLwaaK8VzGzW6iprz2fDODDhG3kAozw+2MAP04mc+wvMPMbetr8UTArt/2ubdsEwO/bYx9PRN8MYBfArTB5xn4JRsD/JBH9AoBfGHgeCQlbRxIWCV1Y2P81mj7StRbAKlQIqc6Z+N2XlIwAfC4zv33NurrKIwBvYeaPDnauXgN9f1V7LDX0z5n5la1KzRomnwrgfxDRS5j5/+tp26uY+Quja2cAfhhGa3s3EX0bmuf2qTCC8TMA/BsiesqK+0hI2BqSzeIGBDNfAvCI4Li/GMBrBi4BgN8B8FwAIKIPBvChHedcgVk32uGvAXwwmZTQFwB8gt3/NgB3EdGT7bYcMF8J4J8Lnv5pI27pt2BST2fWJvHxdv/bAdxBRB9tyyqI6CnM/AiAK0T0Ufa8Lxgou689rwTw1VaDABF9gLVv/C0A9zPzj8GkGn96T7l/AOBjiehv2+t3iegD0AiGB60N4/PscQXg8cz8GwC+GYYSO4f2M09IOBYkzeLGwC4RvUdsfx/M+gP/hYh2AbwTwJetKOOHAbzU0k8uLfQleQIzP0REv2uNxb/KzN9ERD9rz/0Lex2YeU5mydn/TUQPwggit/LaiwD8AIA32QH6HoQUURdeDuAfwKSl/nNYwcfMS0sT/aAVVrkt+y0wKdJ/jIj2Yew0l9rFDrbnv8JQdK+3+x+AWYP6WQC+iYhKAFcB/NOuQpn5ATLG/pcR0dTu/lZm/nMi+jF7L/fA0FqAodB+wt4HAfh+Zr5oDfH/yxre/zkz//aKZ5WQsBFSivKEUSCiDEBhB/onw6xf8AHMvDyFtjwLxqi7SogMlXGOma/a38+HWfPg67bTwtOFFULPYOavPe22JJwdJM0iYSx2AfyGpV0IwFefhqCwWAL4ECL6Fd481uJTiegFMN/AX8N4EF33IKJvgHGn/bnTbkvC2ULSLBISjhlE9Icw7rUSX8zMf3oa7UlI2ARJWCQkJCQkrETyhkpISEhIWIkkLBISEhISViIJi4SEhISElUjCIiEhISFhJf7/5gyz8PumUgwAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ds_avg.tas.plot(label=\"weighted\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Calculate grouped averages\n", "\n", "Related API: [xarray.Dataset.temporal.group_average()](../generated/xarray.Dataset.temporal.group_average.rst)\n", "\n", "Helpful knowledge:\n", "\n", "* Each specified frequency has predefined groups for grouping time coordinates.\n", "\n", " * The table below maps type of averages with its API frequency and grouping convention.\n", "\n", " | Type of Averages | API Frequency | Group By |\n", " |---\t|---\t|---\t|\n", " | Yearly | ``freq=“year”`` | year \t|\n", " | Monthly | ``freq=“month”`` \t| year, month |\n", " | Seasonal | ``freq=“season”`` | year, season |\n", " | Custom seasonal \t| ``freq=\"season\"`` and
``season_config={\"custom_seasons\": <2D ARRAY>}`` | year, season \t|\n", " | Daily | ``freq=“day”`` | year, month, day |\n", " | Hourly | ``freq=“hour”`` | year, month, day, hour |\n", "\n", " * The grouping conventions are based on [CDAT/cdutil](https://github.com/CDAT/cdutil/blob/b823b69db46bb76536db7d435e72075fc3975c65/cdutil/times.py#L1620-L1640), except for daily and hourly means which aren't implemented in CDAT/cdutil.\n", "\n", "* Masked (missing) data is automatically handled.\n", " * The weight of masked (missing) data are excluded when averages are calculated. This is the same as giving them a weight of 0." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Open the ``Dataset``\n", "\n", "In this example, we will be calculating the weighted grouped time averages for `tas` data.\n", "\n", "We are using xarray's OPeNDAP support to read a netCDF4 dataset file directly from its source. The data is not loaded over the network until we perform operations on it (e.g., temperature unit adjustment).\n", "\n", "*More information on the xarray's OPeNDAP support can be found [here](https://docs.xarray.dev/en/stable/user-guide/io.html#opendap).*" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
       "Dimensions:    (time: 1980, bnds: 2, lat: 145, lon: 192)\n",
       "Coordinates:\n",
       "  * time       (time) datetime64[ns] 1850-01-16T12:00:00 ... 2014-12-16T12:00:00\n",
       "  * lat        (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n",
       "  * lon        (lon) float64 0.0 1.875 3.75 5.625 ... 352.5 354.4 356.2 358.1\n",
       "    height     float64 2.0\n",
       "Dimensions without coordinates: bnds\n",
       "Data variables:\n",
       "    time_bnds  (time, bnds) datetime64[ns] 1850-01-01 1850-02-01 ... 2015-01-01\n",
       "    lat_bnds   (lat, bnds) float64 -90.0 -89.38 -89.38 ... 89.38 89.38 90.0\n",
       "    lon_bnds   (lon, bnds) float64 -0.9375 0.9375 0.9375 ... 357.2 357.2 359.1\n",
       "    tas        (time, lat, lon) float32 -27.19 -27.19 -27.19 ... -25.29 -25.29\n",
       "Attributes: (12/49)\n",
       "    Conventions:                     CF-1.7 CMIP-6.2\n",
       "    activity_id:                     CMIP\n",
       "    branch_method:                   standard\n",
       "    branch_time_in_child:            0.0\n",
       "    branch_time_in_parent:           87658.0\n",
       "    creation_date:                   2020-06-05T04:06:11Z\n",
       "    ...                              ...\n",
       "    version:                         v20200605\n",
       "    license:                         CMIP6 model data produced by CSIRO is li...\n",
       "    cmor_version:                    3.4.0\n",
       "    _NCProperties:                   version=2,netcdf=4.6.2,hdf5=1.10.5\n",
       "    tracking_id:                     hdl:21.14100/af78ae5e-f3a6-4e99-8cfe-5f2...\n",
       "    DODS_EXTRA.Unlimited_Dimension:  time
" ], "text/plain": [ "\n", "Dimensions: (time: 1980, bnds: 2, lat: 145, lon: 192)\n", "Coordinates:\n", " * time (time) datetime64[ns] 1850-01-16T12:00:00 ... 2014-12-16T12:00:00\n", " * lat (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n", " * lon (lon) float64 0.0 1.875 3.75 5.625 ... 352.5 354.4 356.2 358.1\n", " height float64 2.0\n", "Dimensions without coordinates: bnds\n", "Data variables:\n", " time_bnds (time, bnds) datetime64[ns] ...\n", " lat_bnds (lat, bnds) float64 ...\n", " lon_bnds (lon, bnds) float64 ...\n", " tas (time, lat, lon) float32 -27.19 -27.19 -27.19 ... -25.29 -25.29\n", "Attributes: (12/49)\n", " Conventions: CF-1.7 CMIP-6.2\n", " activity_id: CMIP\n", " branch_method: standard\n", " branch_time_in_child: 0.0\n", " branch_time_in_parent: 87658.0\n", " creation_date: 2020-06-05T04:06:11Z\n", " ... ...\n", " version: v20200605\n", " license: CMIP6 model data produced by CSIRO is li...\n", " cmor_version: 3.4.0\n", " _NCProperties: version=2,netcdf=4.6.2,hdf5=1.10.5\n", " tracking_id: hdl:21.14100/af78ae5e-f3a6-4e99-8cfe-5f2...\n", " DODS_EXTRA.Unlimited_Dimension: time" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "filepath = \"http://esgf.nci.org.au/thredds/dodsC/master/CMIP6/CMIP/CSIRO/ACCESS-ESM1-5/historical/r10i1p1f1/Amon/tas/gn/v20200605/tas_Amon_ACCESS-ESM1-5_historical_r10i1p1f1_gn_185001-201412.nc\"\n", "\n", "ds = xcdat.open_dataset(filepath)\n", "\n", "# Unit adjust (-273.15, K to C)\n", "ds[\"tas\"] = ds.tas - 273.15\n", "\n", "ds\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Yearly Averages\n", "\n", "**Group time coordinates by year**\n" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "ds_yearly = ds.temporal.group_average(\"tas\", freq=\"year\", weighted=True)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'tas' (time: 165, lat: 145, lon: 192)>\n",
       "array([[[-48.755733, -48.755733, -48.755733, ..., -48.755733,\n",
       "         -48.755733, -48.755733],\n",
       "        [-45.652065, -45.693024, -45.73506 , ..., -45.52128 ,\n",
       "         -45.563866, -45.60669 ],\n",
       "        [-44.775234, -44.905838, -45.03297 , ..., -44.37118 ,\n",
       "         -44.50631 , -44.640503],\n",
       "        ...,\n",
       "        [-20.505976, -20.481321, -20.454565, ..., -20.588959,\n",
       "         -20.557522, -20.530872],\n",
       "        [-20.797592, -20.784252, -20.775455, ..., -20.83268 ,\n",
       "         -20.823357, -20.807684],\n",
       "        [-21.201149, -21.201149, -21.201149, ..., -21.201149,\n",
       "         -21.201149, -21.201149]],\n",
       "\n",
       "       [[-48.95255 , -48.95255 , -48.95255 , ..., -48.95255 ,\n",
       "         -48.95255 , -48.95255 ],\n",
       "        [-45.83191 , -45.864902, -45.89875 , ..., -45.73217 ,\n",
       "         -45.76544 , -45.798595],\n",
       "        [-44.935368, -45.037956, -45.13801 , ..., -44.61143 ,\n",
       "         -44.71986 , -44.829372],\n",
       "...\n",
       "        [-14.916271, -14.899261, -14.88381 , ..., -14.99543 ,\n",
       "         -14.965137, -14.938532],\n",
       "        [-15.405922, -15.396681, -15.385955, ..., -15.432463,\n",
       "         -15.426056, -15.413568],\n",
       "        [-15.945   , -15.945   , -15.945   , ..., -15.945   ,\n",
       "         -15.945   , -15.945   ]],\n",
       "\n",
       "       [[-47.59732 , -47.59732 , -47.59732 , ..., -47.59732 ,\n",
       "         -47.59732 , -47.59732 ],\n",
       "        [-44.721367, -44.763428, -44.803505, ..., -44.592392,\n",
       "         -44.634445, -44.678226],\n",
       "        [-43.85032 , -43.969563, -44.08714 , ..., -43.4709  ,\n",
       "         -43.596764, -43.72408 ],\n",
       "        ...,\n",
       "        [-14.52023 , -14.474079, -14.432307, ..., -14.675514,\n",
       "         -14.620932, -14.567368],\n",
       "        [-14.911236, -14.892309, -14.869016, ..., -14.982012,\n",
       "         -14.962668, -14.938723],\n",
       "        [-15.618406, -15.618406, -15.618406, ..., -15.618406,\n",
       "         -15.618406, -15.618406]]], dtype=float32)\n",
       "Coordinates:\n",
       "  * lat      (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n",
       "  * lon      (lon) float64 0.0 1.875 3.75 5.625 7.5 ... 352.5 354.4 356.2 358.1\n",
       "    height   float64 2.0\n",
       "  * time     (time) object 1850-01-01 00:00:00 ... 2014-01-01 00:00:00\n",
       "Attributes:\n",
       "    operation:  temporal_avg\n",
       "    mode:       group_average\n",
       "    freq:       year\n",
       "    weighted:   True
" ], "text/plain": [ "\n", "array([[[-48.755733, -48.755733, -48.755733, ..., -48.755733,\n", " -48.755733, -48.755733],\n", " [-45.652065, -45.693024, -45.73506 , ..., -45.52128 ,\n", " -45.563866, -45.60669 ],\n", " [-44.775234, -44.905838, -45.03297 , ..., -44.37118 ,\n", " -44.50631 , -44.640503],\n", " ...,\n", " [-20.505976, -20.481321, -20.454565, ..., -20.588959,\n", " -20.557522, -20.530872],\n", " [-20.797592, -20.784252, -20.775455, ..., -20.83268 ,\n", " -20.823357, -20.807684],\n", " [-21.201149, -21.201149, -21.201149, ..., -21.201149,\n", " -21.201149, -21.201149]],\n", "\n", " [[-48.95255 , -48.95255 , -48.95255 , ..., -48.95255 ,\n", " -48.95255 , -48.95255 ],\n", " [-45.83191 , -45.864902, -45.89875 , ..., -45.73217 ,\n", " -45.76544 , -45.798595],\n", " [-44.935368, -45.037956, -45.13801 , ..., -44.61143 ,\n", " -44.71986 , -44.829372],\n", "...\n", " [-14.916271, -14.899261, -14.88381 , ..., -14.99543 ,\n", " -14.965137, -14.938532],\n", " [-15.405922, -15.396681, -15.385955, ..., -15.432463,\n", " -15.426056, -15.413568],\n", " [-15.945 , -15.945 , -15.945 , ..., -15.945 ,\n", " -15.945 , -15.945 ]],\n", "\n", " [[-47.59732 , -47.59732 , -47.59732 , ..., -47.59732 ,\n", " -47.59732 , -47.59732 ],\n", " [-44.721367, -44.763428, -44.803505, ..., -44.592392,\n", " -44.634445, -44.678226],\n", " [-43.85032 , -43.969563, -44.08714 , ..., -43.4709 ,\n", " -43.596764, -43.72408 ],\n", " ...,\n", " [-14.52023 , -14.474079, -14.432307, ..., -14.675514,\n", " -14.620932, -14.567368],\n", " [-14.911236, -14.892309, -14.869016, ..., -14.982012,\n", " -14.962668, -14.938723],\n", " [-15.618406, -15.618406, -15.618406, ..., -15.618406,\n", " -15.618406, -15.618406]]], dtype=float32)\n", "Coordinates:\n", " * lat (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n", " * lon (lon) float64 0.0 1.875 3.75 5.625 7.5 ... 352.5 354.4 356.2 358.1\n", " height float64 2.0\n", " * time (time) object 1850-01-01 00:00:00 ... 2014-01-01 00:00:00\n", "Attributes:\n", " operation: temporal_avg\n", " mode: group_average\n", " freq: year\n", " weighted: True" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds_yearly.tas" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "![tas yearly averages](../examples/temporal-average-yearly.gif)\n", "\n", "*This GIF was created using [xmovie](https://github.com/jbusecke/xmovie).*\n", "\n", "Sample ``xmovie`` code:\n", "```python\n", "import xmovie\n", "mov = xmovie.Movie(ds_yearly_avg.tas)\n", "mov.save(\"temporal-average-yearly.gif\")\n", "```" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Seasonal Averages\n", "**Group time coordinates by year and season**" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "ds_season = ds.temporal.group_average(\"tas\", freq=\"season\", weighted=True)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'tas' (time: 661, lat: 145, lon: 192)>\n",
       "array([[[-32.705883  , -32.705883  , -32.705883  , ..., -32.705883  ,\n",
       "         -32.705883  , -32.705883  ],\n",
       "        [-30.993767  , -31.037586  , -31.089327  , ..., -30.845623  ,\n",
       "         -30.894127  , -30.94401   ],\n",
       "        [-30.02515   , -30.145437  , -30.26419   , ..., -29.660372  ,\n",
       "         -29.78108   , -29.902878  ],\n",
       "        ...,\n",
       "        [-37.72314   , -37.685493  , -37.654167  , ..., -37.8262    ,\n",
       "         -37.790344  , -37.75683   ],\n",
       "        [-38.274647  , -38.263725  , -38.250145  , ..., -38.292183  ,\n",
       "         -38.290638  , -38.28456   ],\n",
       "        [-38.743587  , -38.743587  , -38.743587  , ..., -38.743587  ,\n",
       "         -38.743587  , -38.743587  ]],\n",
       "\n",
       "       [[-54.290863  , -54.290863  , -54.290863  , ..., -54.290863  ,\n",
       "         -54.290863  , -54.290863  ],\n",
       "        [-51.117714  , -51.175236  , -51.230553  , ..., -50.935165  ,\n",
       "         -50.99657   , -51.056145  ],\n",
       "        [-50.318047  , -50.486664  , -50.649567  , ..., -49.79003   ,\n",
       "         -49.970078  , -50.14521   ],\n",
       "...\n",
       "        [-12.342774  , -12.2246685 , -12.106632  , ..., -12.744922  ,\n",
       "         -12.609088  , -12.478392  ],\n",
       "        [-13.126404  , -13.066109  , -13.003876  , ..., -13.306077  ,\n",
       "         -13.258715  , -13.19972   ],\n",
       "        [-14.288469  , -14.288469  , -14.288469  , ..., -14.288469  ,\n",
       "         -14.288469  , -14.288469  ]],\n",
       "\n",
       "       [[-28.990494  , -28.990494  , -28.990494  , ..., -28.990494  ,\n",
       "         -28.990494  , -28.990494  ],\n",
       "        [-28.192917  , -28.224579  , -28.261307  , ..., -28.095932  ,\n",
       "         -28.125992  , -28.15802   ],\n",
       "        [-27.607407  , -27.705643  , -27.805115  , ..., -27.311615  ,\n",
       "         -27.410828  , -27.508362  ],\n",
       "        ...,\n",
       "        [-24.256271  , -24.140594  , -24.037537  , ..., -24.61853   ,\n",
       "         -24.488495  , -24.36644   ],\n",
       "        [-24.629013  , -24.613388  , -24.549866  , ..., -24.752045  ,\n",
       "         -24.721603  , -24.666412  ],\n",
       "        [-25.28923   , -25.28923   , -25.28923   , ..., -25.28923   ,\n",
       "         -25.28923   , -25.28923   ]]], dtype=float32)\n",
       "Coordinates:\n",
       "  * lat      (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n",
       "  * lon      (lon) float64 0.0 1.875 3.75 5.625 7.5 ... 352.5 354.4 356.2 358.1\n",
       "    height   float64 2.0\n",
       "  * time     (time) object 1850-01-01 00:00:00 ... 2015-01-01 00:00:00\n",
       "Attributes:\n",
       "    operation:            temporal_avg\n",
       "    mode:                 group_average\n",
       "    freq:                 season\n",
       "    weighted:             True\n",
       "    dec_mode:             DJF\n",
       "    drop_incomplete_djf:  False
" ], "text/plain": [ "\n", "array([[[-32.705883 , -32.705883 , -32.705883 , ..., -32.705883 ,\n", " -32.705883 , -32.705883 ],\n", " [-30.993767 , -31.037586 , -31.089327 , ..., -30.845623 ,\n", " -30.894127 , -30.94401 ],\n", " [-30.02515 , -30.145437 , -30.26419 , ..., -29.660372 ,\n", " -29.78108 , -29.902878 ],\n", " ...,\n", " [-37.72314 , -37.685493 , -37.654167 , ..., -37.8262 ,\n", " -37.790344 , -37.75683 ],\n", " [-38.274647 , -38.263725 , -38.250145 , ..., -38.292183 ,\n", " -38.290638 , -38.28456 ],\n", " [-38.743587 , -38.743587 , -38.743587 , ..., -38.743587 ,\n", " -38.743587 , -38.743587 ]],\n", "\n", " [[-54.290863 , -54.290863 , -54.290863 , ..., -54.290863 ,\n", " -54.290863 , -54.290863 ],\n", " [-51.117714 , -51.175236 , -51.230553 , ..., -50.935165 ,\n", " -50.99657 , -51.056145 ],\n", " [-50.318047 , -50.486664 , -50.649567 , ..., -49.79003 ,\n", " -49.970078 , -50.14521 ],\n", "...\n", " [-12.342774 , -12.2246685 , -12.106632 , ..., -12.744922 ,\n", " -12.609088 , -12.478392 ],\n", " [-13.126404 , -13.066109 , -13.003876 , ..., -13.306077 ,\n", " -13.258715 , -13.19972 ],\n", " [-14.288469 , -14.288469 , -14.288469 , ..., -14.288469 ,\n", " -14.288469 , -14.288469 ]],\n", "\n", " [[-28.990494 , -28.990494 , -28.990494 , ..., -28.990494 ,\n", " -28.990494 , -28.990494 ],\n", " [-28.192917 , -28.224579 , -28.261307 , ..., -28.095932 ,\n", " -28.125992 , -28.15802 ],\n", " [-27.607407 , -27.705643 , -27.805115 , ..., -27.311615 ,\n", " -27.410828 , -27.508362 ],\n", " ...,\n", " [-24.256271 , -24.140594 , -24.037537 , ..., -24.61853 ,\n", " -24.488495 , -24.36644 ],\n", " [-24.629013 , -24.613388 , -24.549866 , ..., -24.752045 ,\n", " -24.721603 , -24.666412 ],\n", " [-25.28923 , -25.28923 , -25.28923 , ..., -25.28923 ,\n", " -25.28923 , -25.28923 ]]], dtype=float32)\n", "Coordinates:\n", " * lat (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n", " * lon (lon) float64 0.0 1.875 3.75 5.625 7.5 ... 352.5 354.4 356.2 358.1\n", " height float64 2.0\n", " * time (time) object 1850-01-01 00:00:00 ... 2015-01-01 00:00:00\n", "Attributes:\n", " operation: temporal_avg\n", " mode: group_average\n", " freq: season\n", " weighted: True\n", " dec_mode: DJF\n", " drop_incomplete_djf: False" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds_season.tas" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Notice that the season of each time coordinate is represented by its middle month.**\n", "\n", "- \"DJF\" is represented by month 1 (\"J\"/January)\n", "- \"MAM\" is represented by month 4 (\"A\"/April)\n", "- \"JJA\" is represented by month 7 (\"J\"/July)\n", "- \"SON\" is represented by month 10 (\"O\"/October).\n", "\n", "This is implementation design was used because `datetime` objects do not distinguish seasons, so the middle month is used instead.\n" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'time' (time: 661)>\n",
       "array([cftime.DatetimeProlepticGregorian(1850, 1, 1, 0, 0, 0, 0, has_year_zero=True),\n",
       "       cftime.DatetimeProlepticGregorian(1850, 4, 1, 0, 0, 0, 0, has_year_zero=True),\n",
       "       cftime.DatetimeProlepticGregorian(1850, 7, 1, 0, 0, 0, 0, has_year_zero=True),\n",
       "       ...,\n",
       "       cftime.DatetimeProlepticGregorian(2014, 7, 1, 0, 0, 0, 0, has_year_zero=True),\n",
       "       cftime.DatetimeProlepticGregorian(2014, 10, 1, 0, 0, 0, 0, has_year_zero=True),\n",
       "       cftime.DatetimeProlepticGregorian(2015, 1, 1, 0, 0, 0, 0, has_year_zero=True)],\n",
       "      dtype=object)\n",
       "Coordinates:\n",
       "    height   float64 2.0\n",
       "  * time     (time) object 1850-01-01 00:00:00 ... 2015-01-01 00:00:00\n",
       "Attributes:\n",
       "    bounds:         time_bnds\n",
       "    axis:           T\n",
       "    long_name:      time\n",
       "    standard_name:  time\n",
       "    _ChunkSizes:    1
" ], "text/plain": [ "\n", "array([cftime.DatetimeProlepticGregorian(1850, 1, 1, 0, 0, 0, 0, has_year_zero=True),\n", " cftime.DatetimeProlepticGregorian(1850, 4, 1, 0, 0, 0, 0, has_year_zero=True),\n", " cftime.DatetimeProlepticGregorian(1850, 7, 1, 0, 0, 0, 0, has_year_zero=True),\n", " ...,\n", " cftime.DatetimeProlepticGregorian(2014, 7, 1, 0, 0, 0, 0, has_year_zero=True),\n", " cftime.DatetimeProlepticGregorian(2014, 10, 1, 0, 0, 0, 0, has_year_zero=True),\n", " cftime.DatetimeProlepticGregorian(2015, 1, 1, 0, 0, 0, 0, has_year_zero=True)],\n", " dtype=object)\n", "Coordinates:\n", " height float64 2.0\n", " * time (time) object 1850-01-01 00:00:00 ... 2015-01-01 00:00:00\n", "Attributes:\n", " bounds: time_bnds\n", " axis: T\n", " long_name: time\n", " standard_name: time\n", " _ChunkSizes: 1" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds_season.time" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Monthly Averages\n", "\n", "**Group time coordinates by year and month**\n", "\n", "For this example, we will be loading a subset of daily time series data for `tas` using OPeNDAP.\n", "\n", "---\n", "\n", "**NOTE:**\n", "\n", "For OPeNDAP servers, the default file size request limit is 500MB in the TDS server configuration. Opening up a dataset over OPeNDAP also introduces an overhead compared to direct file access.\n", "\n", "**The workaround is to use Dask to request the data in manageable chunks, which overcomes file size limitations and can improve performance.**\n", "\n", "We have a few ways to chunk our request:\n", "\n", "1. Specify `chunks` with `\"auto\"` to let Dask determine the chunksize.\n", "2. Specify a specify the file size to chunk on (e.g., `\"100MB\"`) or number of chunks as an integer (`100` for 100 chunks).\n", "\n", "Visit this page to learn more about chunking and performance: https://docs.xarray.dev/en/stable/user-guide/dask.html#chunking-and-performance\n", "\n", "---" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
       "Dimensions:    (time: 18262, bnds: 2, lat: 145, lon: 192)\n",
       "Coordinates:\n",
       "  * time       (time) datetime64[ns] 1850-01-01T12:00:00 ... 1899-12-31T12:00:00\n",
       "  * lat        (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n",
       "  * lon        (lon) float64 0.0 1.875 3.75 5.625 ... 352.5 354.4 356.2 358.1\n",
       "    height     float64 ...\n",
       "Dimensions without coordinates: bnds\n",
       "Data variables:\n",
       "    time_bnds  (time, bnds) datetime64[ns] dask.array<chunksize=(18262, 2), meta=np.ndarray>\n",
       "    lat_bnds   (lat, bnds) float64 dask.array<chunksize=(145, 2), meta=np.ndarray>\n",
       "    lon_bnds   (lon, bnds) float64 dask.array<chunksize=(192, 2), meta=np.ndarray>\n",
       "    tas        (time, lat, lon) float32 dask.array<chunksize=(794, 145, 192), meta=np.ndarray>\n",
       "Attributes: (12/48)\n",
       "    Conventions:                     CF-1.7 CMIP-6.2\n",
       "    activity_id:                     CMIP\n",
       "    branch_method:                   standard\n",
       "    branch_time_in_child:            0.0\n",
       "    branch_time_in_parent:           21915.0\n",
       "    creation_date:                   2019-11-15T17:30:04Z\n",
       "    ...                              ...\n",
       "    variant_label:                   r1i1p1f1\n",
       "    version:                         v20191115\n",
       "    cmor_version:                    3.4.0\n",
       "    tracking_id:                     hdl:21.14100/a9d8ba3a-bcbf-4d54-9970-cfc...\n",
       "    license:                         CMIP6 model data produced by CSIRO is li...\n",
       "    DODS_EXTRA.Unlimited_Dimension:  time
" ], "text/plain": [ "\n", "Dimensions: (time: 18262, bnds: 2, lat: 145, lon: 192)\n", "Coordinates:\n", " * time (time) datetime64[ns] 1850-01-01T12:00:00 ... 1899-12-31T12:00:00\n", " * lat (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n", " * lon (lon) float64 0.0 1.875 3.75 5.625 ... 352.5 354.4 356.2 358.1\n", " height float64 ...\n", "Dimensions without coordinates: bnds\n", "Data variables:\n", " time_bnds (time, bnds) datetime64[ns] dask.array\n", " lat_bnds (lat, bnds) float64 dask.array\n", " lon_bnds (lon, bnds) float64 dask.array\n", " tas (time, lat, lon) float32 dask.array\n", "Attributes: (12/48)\n", " Conventions: CF-1.7 CMIP-6.2\n", " activity_id: CMIP\n", " branch_method: standard\n", " branch_time_in_child: 0.0\n", " branch_time_in_parent: 21915.0\n", " creation_date: 2019-11-15T17:30:04Z\n", " ... ...\n", " variant_label: r1i1p1f1\n", " version: v20191115\n", " cmor_version: 3.4.0\n", " tracking_id: hdl:21.14100/a9d8ba3a-bcbf-4d54-9970-cfc...\n", " license: CMIP6 model data produced by CSIRO is li...\n", " DODS_EXTRA.Unlimited_Dimension: time" ] }, "execution_count": 13, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# The size of this file is approximately 1.45 GB, so we will be chunking our\n", "# request using Dask to avoid hitting the OPeNDAP file size request limit for\n", "# this ESGF node.\n", "ds2 = xcdat.open_dataset(\n", " \"http://esgf-data3.diasjp.net/thredds/dodsC/esg_dataroot/CMIP6/CMIP/CSIRO/ACCESS-ESM1-5/historical/r1i1p1f1/day/tas/gn/v20191115/tas_day_ACCESS-ESM1-5_historical_r1i1p1f1_gn_18500101-18991231.nc\",\n", " chunks={\"time\": \"auto\"},\n", ")\n", "\n", "# Unit adjust (-273.15, K to C)\n", "ds2[\"tas\"] = ds2.tas - 273.15\n", "\n", "ds2" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "ds2_monthly_avg = ds2.temporal.group_average(\"tas\", freq=\"month\", weighted=True)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'tas' (time: 600, lat: 145, lon: 192)>\n",
       "dask.array<stack, shape=(600, 145, 192), dtype=float64, chunksize=(1, 145, 192), chunktype=numpy.ndarray>\n",
       "Coordinates:\n",
       "  * lat      (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n",
       "  * lon      (lon) float64 0.0 1.875 3.75 5.625 7.5 ... 352.5 354.4 356.2 358.1\n",
       "    height   float64 ...\n",
       "  * time     (time) object 1850-01-01 00:00:00 ... 1899-12-01 00:00:00\n",
       "Attributes:\n",
       "    operation:  temporal_avg\n",
       "    mode:       group_average\n",
       "    freq:       month\n",
       "    weighted:   True
" ], "text/plain": [ "\n", "dask.array\n", "Coordinates:\n", " * lat (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n", " * lon (lon) float64 0.0 1.875 3.75 5.625 7.5 ... 352.5 354.4 356.2 358.1\n", " height float64 ...\n", " * time (time) object 1850-01-01 00:00:00 ... 1899-12-01 00:00:00\n", "Attributes:\n", " operation: temporal_avg\n", " mode: group_average\n", " freq: month\n", " weighted: True" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds2_monthly_avg.tas" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Daily Averages\n", "\n", "**Group time coordinates by year, month, and day**\n", "\n", "For this example, we will be opening a subset of 3hr time series data for `tas` using OPeNDAP." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "# The size of this file is approximately 1.17 GB, so we will be chunking our\n", "# request using Dask to avoid hitting the OPeNDAP file size request limit for\n", "# this ESGF node.\n", "ds3 = xcdat.open_dataset(\n", " \"http://esgf.nci.org.au/thredds/dodsC/master/CMIP6/CMIP/CSIRO/ACCESS-ESM1-5/historical/r10i1p1f1/3hr/tas/gn/v20200605/tas_3hr_ACCESS-ESM1-5_historical_r10i1p1f1_gn_201001010300-201501010000.nc\",\n", " chunks={\"time\": \"auto\"}\n", ")\n", "\n", "# Unit adjust (-273.15, K to C)\n", "ds3[\"tas\"] = ds3.tas - 273.15" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'tas' (time: 14608, lat: 145, lon: 192)>\n",
       "dask.array<sub, shape=(14608, 145, 192), dtype=float32, chunksize=(913, 145, 192), chunktype=numpy.ndarray>\n",
       "Coordinates:\n",
       "  * time     (time) datetime64[ns] 2010-01-01T03:00:00 ... 2015-01-01\n",
       "  * lat      (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n",
       "  * lon      (lon) float64 0.0 1.875 3.75 5.625 7.5 ... 352.5 354.4 356.2 358.1\n",
       "    height   float64 ...
" ], "text/plain": [ "\n", "dask.array\n", "Coordinates:\n", " * time (time) datetime64[ns] 2010-01-01T03:00:00 ... 2015-01-01\n", " * lat (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n", " * lon (lon) float64 0.0 1.875 3.75 5.625 7.5 ... 352.5 354.4 356.2 358.1\n", " height float64 ..." ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds3.tas" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "ds3_day_avg = ds3.temporal.group_average(\"tas\", freq=\"day\", weighted=True)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'tas' (time: 1827, lat: 145, lon: 192)>\n",
       "dask.array<stack, shape=(1827, 145, 192), dtype=float64, chunksize=(1, 145, 192), chunktype=numpy.ndarray>\n",
       "Coordinates:\n",
       "  * lat      (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n",
       "  * lon      (lon) float64 0.0 1.875 3.75 5.625 7.5 ... 352.5 354.4 356.2 358.1\n",
       "    height   float64 ...\n",
       "  * time     (time) object 2010-01-01 00:00:00 ... 2015-01-01 00:00:00\n",
       "Attributes:\n",
       "    operation:  temporal_avg\n",
       "    mode:       group_average\n",
       "    freq:       day\n",
       "    weighted:   True
" ], "text/plain": [ "\n", "dask.array\n", "Coordinates:\n", " * lat (lat) float64 -90.0 -88.75 -87.5 -86.25 ... 86.25 87.5 88.75 90.0\n", " * lon (lon) float64 0.0 1.875 3.75 5.625 7.5 ... 352.5 354.4 356.2 358.1\n", " height float64 ...\n", " * time (time) object 2010-01-01 00:00:00 ... 2015-01-01 00:00:00\n", "Attributes:\n", " operation: temporal_avg\n", " mode: group_average\n", " freq: day\n", " weighted: True" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds3_day_avg.tas" ] } ], "metadata": { "anaconda-cloud": {}, "interpreter": { "hash": "937205ea97caa5d37934716e0a0c6b9e4ffcdaf6e0f20ca1ff82361fb532cef2" }, "kernelspec": { "display_name": "Python 3.9.12 ('xcdat_dev')", "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.9.13" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": true, "toc_position": {}, "toc_section_display": true, "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 4 }