From 3617649f97a4d43be53a38ef7462ca1bc5b5086e Mon Sep 17 00:00:00 2001 From: Wolfgang Preimesberger Date: Tue, 22 Aug 2023 10:02:01 +0200 Subject: [PATCH] Metrics for temporal subgroups (#266) * Update environment * First implementation of flexible set metrics * Fix keyword for metrics calculation when reference dataset must be included * Update tests * Update CHANGELOG.rst * Update CHANGELOG.rst * Update env * Remove unnecessary checks for data availability * Undo * Fix Test * Make bootstrapping settings better accessible when using the validation framework * Renamed GenericDatetime to YearlessDatetime and moved to grouping module * Update notebook to include subset metrics and reader adapters * Update tests * Change yearless date name * Fix tests --- docs/examples/validation_framework.ipynb | 546 +++++++++++++++--- src/pytesmo/time_series/grouping.py | 245 +++++++- src/pytesmo/validation_framework/adapters.py | 58 +- .../metric_calculators_adapters.py | 180 +++--- .../validation_framework/validation.py | 7 +- .../test_metric_calculators_adapters.py | 131 +++++ 6 files changed, 980 insertions(+), 187 deletions(-) create mode 100644 tests/test_validation_framework/test_metric_calculators_adapters.py diff --git a/docs/examples/validation_framework.ipynb b/docs/examples/validation_framework.ipynb index aebc5657..5b41390b 100644 --- a/docs/examples/validation_framework.ipynb +++ b/docs/examples/validation_framework.ipynb @@ -74,12 +74,15 @@ "metadata": {}, "outputs": [], "source": [ + "import os\n", + "import numpy as np\n", "from pathlib import Path\n", - "\n", + "from pprint import pprint\n", "import pytesmo.validation_framework.metric_calculators as metrics_calculators\n", "from datetime import datetime\n", "import warnings\n", "from ascat.read_native.cdr import AscatGriddedNcTs\n", + "\n", "from ismn.interface import ISMN_Interface # install ismn: 'pip install ismn'\n", "from pytesmo.validation_framework.validation import Validation\n", "from pytesmo.validation_framework.results_manager import netcdf_results_manager\n", @@ -104,13 +107,13 @@ "name": "stdout", "output_type": "stream", "text": [ - "Data is stored in: /tmp/tmprnczvqbv\n" + "Data is stored in: /tmp/tmppeisd__e\n" ] } ], "source": [ "from tempfile import mkdtemp\n", - "output_folder = mkdtemp()\n", + "output_folder = Path(mkdtemp())\n", "print('Data is stored in:', output_folder)" ] }, @@ -146,6 +149,7 @@ " grid_filename=ascat_grid_fname,\n", " static_layer_path=static_layer_path\n", " )\n", + "\n", "ascat_reader.read_bulk = True" ] }, @@ -165,7 +169,7 @@ "name": "stdout", "output_type": "stream", "text": [ - "Processing metadata for all ismn stations into folder /home/samuel/projects/QA4SM/pytesmo/tests/test-data/ismn/multinetwork/header_values.\n", + "Processing metadata for all ismn stations into folder /home/wpreimes/shares/home/code/pytesmo/tests/test-data/ismn/multinetwork/header_values.\n", "This may take a few minutes, but is only done once...\n", "Hint: Use `parallel=True` to speed up metadata generation for large datasets\n" ] @@ -174,7 +178,7 @@ "name": "stderr", "output_type": "stream", "text": [ - "Files Processed: 100%|█████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████████| 8/8 [00:00<00:00, 226.95it/s]" + "Files Processed: 100%|██████████| 8/8 [00:00<00:00, 65.44it/s]" ] }, { @@ -182,8 +186,8 @@ "output_type": "stream", "text": [ "Metadata generation finished after 0 Seconds.\n", - "Metadata and Log stored in /tmp/tmpyw8mctbb\n", - "Found existing ismn metadata in /tmp/tmpyw8mctbb/header_values.csv.\n" + "Metadata and Log stored in /tmp/tmp_72gspsz\n", + "Found existing ismn metadata in /tmp/tmp_72gspsz/header_values.csv.\n" ] }, { @@ -212,11 +216,6 @@ "**DO NOT CHANGE** the name ***jobs*** because it will be searched during the parallel processing!" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [] - }, { "cell_type": "code", "execution_count": 5, @@ -332,26 +331,26 @@ }, { "cell_type": "code", - "execution_count": 8, + "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ - "/home/samuel/projects/QA4SM/pytesmo/src/pytesmo/validation_framework/validation.py:141: UserWarning: You are using the default temporal matcher. If you are using one of the newer metric calculators (PairwiseIntercomparisonMetrics, TripleCollocationMetrics) you should probably use `make_combined_temporal_matcher` instead. Have a look at the documentation of the metric calculators for more info.\n", + "/home/wpreimes/shares/home/code/pytesmo/src/pytesmo/validation_framework/validation.py:144: UserWarning: You are using the default temporal matcher. If you are using one of the newer metric calculators (PairwiseIntercomparisonMetrics, TripleCollocationMetrics) you should probably use `make_combined_temporal_matcher` instead. Have a look at the documentation of the metric calculators for more info.\n", " warnings.warn(\n" ] } ], "source": [ "period = [datetime(2007, 1, 1), datetime(2014, 12, 31)]\n", - "basic_metrics = metrics_calculators.BasicMetrics(other_name='k1')\n", + "basic_metrics = metrics_calculators.PairwiseIntercomparisonMetrics()\n", "\n", "process = Validation(\n", " datasets, 'ISMN',\n", " temporal_ref='ASCAT',\n", - " scaling='cdf_match',\n", + " scaling='mean_std',\n", " scaling_ref='ASCAT',\n", " metrics_calculators={(2, 2): basic_metrics.calc_metrics},\n", " period=period)" @@ -375,7 +374,7 @@ "together and then combinations of two input datasets are given to one metric calculator while all three datasets\n", "are given to another metric calculator. This could look like this:\n", "\n", - "```python\n", + "```\n", "{ (3 ,2): metric_calc,\n", " (3, 3): triple_collocation}\n", "```\n", @@ -386,44 +385,52 @@ }, { "cell_type": "code", - "execution_count": 9, + "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "{(('ASCAT', 'sm'), ('ISMN', 'soil_moisture')): {'BIAS': array([0.18759985], dtype=float32),\n", - " 'R': array([0.47605154], dtype=float32),\n", - " 'RMSD': array([11.537956], dtype=float32),\n", + "{(('ASCAT', 'sm'), ('ISMN', 'soil_moisture')): {'BIAS': array([1.4210855e-14], dtype=float32),\n", + " 'BIAS_ci_lower': array([-1.124606], dtype=float32),\n", + " 'BIAS_ci_upper': array([1.124606], dtype=float32),\n", + " 'R': array([0.5342869], dtype=float32),\n", + " 'RMSD': array([10.789426], dtype=float32),\n", + " 'RSS': array([41558.98], dtype=float32),\n", + " 'R_ci_lower': array([0.45576647], dtype=float32),\n", + " 'R_ci_upper': array([0.60455596], dtype=float32),\n", " 'gpi': array([0], dtype=int32),\n", " 'lat': array([33.8833]),\n", " 'lon': array([102.1333]),\n", + " 'mse': array([116.41171], dtype=float32),\n", + " 'mse_bias': array([2.019484e-28], dtype=float32),\n", + " 'mse_corr': array([116.41171], dtype=float32),\n", + " 'mse_var': array([3.1554436e-30], dtype=float32),\n", " 'n_obs': array([357], dtype=int32),\n", - " 'p_R': array([1.3616014e-21], dtype=float32),\n", + " 'p_R': array([9.664754e-28], dtype=float32),\n", " 'p_rho': array([2.471651e-28], dtype=float32),\n", - " 'p_tau': array([nan], dtype=float32),\n", + " 'p_tau': array([2.1020434e-26], dtype=float32),\n", " 'rho': array([0.53934574], dtype=float32),\n", - " 'tau': array([nan], dtype=float32)}}\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/home/samuel/projects/QA4SM/pytesmo/src/pytesmo/cdf_matching.py:174: UserWarning: The bins have been resized\n", - " warnings.warn(\"The bins have been resized\")\n" + " 'rho_ci_lower': array([0.45559874], dtype=float32),\n", + " 'rho_ci_upper': array([0.61362934], dtype=float32),\n", + " 'status': array([0], dtype=int32),\n", + " 'tau': array([0.3907124], dtype=float32),\n", + " 'tau_ci_lower': array([0.35201582], dtype=float32),\n", + " 'tau_ci_upper': array([0.42807564], dtype=float32),\n", + " 'urmsd': array([10.789426], dtype=float32),\n", + " 'urmsd_ci_lower': array([10.065898], dtype=float32),\n", + " 'urmsd_ci_upper': array([11.661137], dtype=float32)}}\n" ] } ], "source": [ - "save_path = output_folder\n", + "save_path = output_folder / 'ascat_ismn'\n", "\n", - "import pprint\n", "for job in jobs:\n", - " \n", + "\n", " results = process.calc(*job)\n", - " pprint.pprint(results)\n", + " pprint(results)\n", " netcdf_results_manager(results, save_path)" ] }, @@ -444,7 +451,7 @@ }, { "cell_type": "code", - "execution_count": 10, + "execution_count": 12, "metadata": {}, "outputs": [ { @@ -457,15 +464,32 @@ "_row_size [0]\n", "time []\n", "gpi [0]\n", - "n_obs [357]\n", - "R [0.47605154]\n", - "p_R [1.3616014e-21]\n", + "status [0]\n", + "R [0.5342869]\n", + "p_R [9.664754e-28]\n", + "BIAS [1.4210855e-14]\n", + "RMSD [10.789426]\n", + "mse [116.41171]\n", + "RSS [41558.98]\n", + "mse_corr [116.41171]\n", + "mse_bias [2.019484e-28]\n", + "urmsd [10.789426]\n", + "mse_var [3.1554436e-30]\n", "rho [0.53934574]\n", "p_rho [2.471651e-28]\n", - "RMSD [11.537956]\n", - "BIAS [0.18759985]\n", - "tau [nan]\n", - "p_tau [nan]\n" + "tau [0.3907124]\n", + "p_tau [2.1020434e-26]\n", + "BIAS_ci_lower [-1.124606]\n", + "BIAS_ci_upper [1.124606]\n", + "urmsd_ci_lower [10.065898]\n", + "urmsd_ci_upper [11.661137]\n", + "R_ci_lower [0.45576647]\n", + "R_ci_upper [0.60455596]\n", + "rho_ci_lower [0.45559874]\n", + "rho_ci_upper [0.61362934]\n", + "tau_ci_lower [0.35201582]\n", + "tau_ci_upper [0.42807564]\n", + "n_obs [357]\n" ] } ], @@ -478,6 +502,174 @@ " print(var, ds.variables[var][:])" ] }, + { + "cell_type": "markdown", + "metadata": { + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "## Metric Calculator Adapters\n", + "\n", + "Metric calculators compute a set of comparison metrics based on the passed input. Usually the input are two or more collocated time series. However, to validate soil moisture it is often desired to assess the quality for certain temporal subsets independently. E.g. if a varying level of agreement between two datasets is expected for two different seasons (e.g. wet/dry season, summer/winter) or if the varying level of agreement should be assessed (e.g. on a seasonal/monthly basis). In this case pytesmo provides adapters to split up the input time series before computing validation metrics. These adapters work with any (properly implemented) metric calculator. You can use one of the predefined adapters in `pytesmo.validation_framework.metric_calculators_adapters`. Currently, there are 2 options:\n", + "\n", + "- `SubsetsMetricsAdapter`: This adapter lets you define arbitrary temporal subsets of your time series data before metrics computation. You can compute metrics for a certain set of datetimes or datetime ranges.\n", + " - `MonthsMetricsAdapter`: This is a version of a SubsetsMetricsAdapter that allows splitting up the data based on their month. You can e.g. compute validation metrics for each month of a year individually, or for each season (i.e. 3 * 4 months), or any other combination of months.\n", + "\n", + "Let's look at an example. Here we use the base `SubsetsMetricsAdapter`. We compute the same metrics as before, but apply the metrics calculator to a set of predefined temporal subgroups. Note that groups in this example were chosen to demonstrate the feature range of the `TsDistributor` rather than to produce meaningful results.\n", + "\n", + "We create a SubsetsMetricsAdapter that takes the original metrics calculator and a list of named, temporal subsets (dict).\n" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": { + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [], + "source": [ + "import os\n", + "from pytesmo.validation_framework.metric_calculators_adapters import SubsetsMetricsAdapter, TsDistributor\n", + "from pytesmo.time_series.grouping import YearlessDatetime\n", + "from datetime import datetime\n", + "import shutil\n", + "\n", + "subsets = {\n", + " # The first subset does only include values in August 2010 and 2011 + July 31st of both years\n", + " 'August2010/11': TsDistributor(dates=[datetime(year=2009, month=7, day=31), datetime(2010, 7, 31)],\n", + " date_ranges=[(datetime(2010, 7, 1, 0, 0), datetime(2010, 7, 31)),\n", + " (datetime(2011, 7, 1), datetime(2011, 7, 31))]),\n", + " \n", + " # The second subset includes all values from June to September of ANY YEAR, but not August\n", + " 'JJS': TsDistributor(yearless_date_ranges=[\n", + " (YearlessDatetime(month=6, day=1, hour=0, minute=0), YearlessDatetime(7, 31)),\n", + " (YearlessDatetime(9, 1), YearlessDatetime(9, 30))\n", + " ]),\n", + " \n", + " # The third group includes all values from Feb 28th to the end of April of ANY YEAR\n", + " 'MarchApril': TsDistributor(yearless_date_ranges=[(YearlessDatetime(2, 28), YearlessDatetime(4, 30))]),\n", + "}\n", + "\n", + "adapted_intercomparison_metrics = SubsetsMetricsAdapter(\n", + " metrics_calculators.PairwiseIntercomparisonMetrics(calc_kendall=False, calc_spearman=False),\n", + " subsets,\n", + " group_results='join')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "The adapted metrics calculator can now be used as before. The results will contain the chosen group names accordingly, as defined in the adapter. The results manager will write them to a netcdf file. Note that this only works as `group_results='join'` was selected, otherwise the output format would different from what the results manager expects." + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": { + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{(('ASCAT', 'sm'), ('ISMN', 'soil_moisture')): {'August2010/11|BIAS': array([0.02893662], dtype=float32),\n", + " 'August2010/11|BIAS_ci_lower': array([0.01474692], dtype=float32),\n", + " 'August2010/11|BIAS_ci_upper': array([0.04312633], dtype=float32),\n", + " 'August2010/11|R': array([0.7176565], dtype=float32),\n", + " 'August2010/11|RMSD': array([0.04613708], dtype=float32),\n", + " 'August2010/11|RSS': array([0.05960163], dtype=float32),\n", + " 'August2010/11|R_ci_lower': array([0.47057065], dtype=float32),\n", + " 'August2010/11|R_ci_upper': array([0.8603755], dtype=float32),\n", + " 'August2010/11|mse': array([0.00212863], dtype=float32),\n", + " 'August2010/11|mse_bias': array([0.00083733], dtype=float32),\n", + " 'August2010/11|mse_corr': array([0.00128012], dtype=float32),\n", + " 'August2010/11|mse_var': array([1.1178111e-05], dtype=float32),\n", + " 'August2010/11|n_obs': array([28], dtype=int32),\n", + " 'August2010/11|p_R': array([1.718043e-05], dtype=float32),\n", + " 'August2010/11|status': array([0], dtype=int32),\n", + " 'August2010/11|urmsd': array([0.03593468], dtype=float32),\n", + " 'August2010/11|urmsd_ci_lower': array([0.02893201], dtype=float32),\n", + " 'August2010/11|urmsd_ci_upper': array([0.04980955], dtype=float32),\n", + " 'JJS|BIAS': array([0.02439872], dtype=float32),\n", + " 'JJS|BIAS_ci_lower': array([0.01300029], dtype=float32),\n", + " 'JJS|BIAS_ci_upper': array([0.03579714], dtype=float32),\n", + " 'JJS|R': array([0.43224052], dtype=float32),\n", + " 'JJS|RMSD': array([0.08166534], dtype=float32),\n", + " 'JJS|RSS': array([1.2204688], dtype=float32),\n", + " 'JJS|R_ci_lower': array([0.3063946], dtype=float32),\n", + " 'JJS|R_ci_upper': array([0.54323655], dtype=float32),\n", + " 'JJS|mse': array([0.00666923], dtype=float32),\n", + " 'JJS|mse_bias': array([0.0005953], dtype=float32),\n", + " 'JJS|mse_corr': array([0.00590448], dtype=float32),\n", + " 'JJS|mse_var': array([0.00016945], dtype=float32),\n", + " 'JJS|n_obs': array([183], dtype=int32),\n", + " 'JJS|p_R': array([9.96128e-10], dtype=float32),\n", + " 'JJS|status': array([0], dtype=int32),\n", + " 'JJS|urmsd': array([0.07793543], dtype=float32),\n", + " 'JJS|urmsd_ci_lower': array([0.07087912], dtype=float32),\n", + " 'JJS|urmsd_ci_upper': array([0.0870943], dtype=float32),\n", + " 'MarchApril|BIAS': array([-0.0636023], dtype=float32),\n", + " 'MarchApril|BIAS_ci_lower': array([-0.07745089], dtype=float32),\n", + " 'MarchApril|BIAS_ci_upper': array([-0.0497537], dtype=float32),\n", + " 'MarchApril|R': array([0.8920552], dtype=float32),\n", + " 'MarchApril|RMSD': array([0.07466082], dtype=float32),\n", + " 'MarchApril|RSS': array([0.1895241], dtype=float32),\n", + " 'MarchApril|R_ci_lower': array([0.7931545], dtype=float32),\n", + " 'MarchApril|R_ci_upper': array([0.94511515], dtype=float32),\n", + " 'MarchApril|mse': array([0.00557424], dtype=float32),\n", + " 'MarchApril|mse_bias': array([0.00404525], dtype=float32),\n", + " 'MarchApril|mse_corr': array([0.00138048], dtype=float32),\n", + " 'MarchApril|mse_var': array([0.0001485], dtype=float32),\n", + " 'MarchApril|n_obs': array([34], dtype=int32),\n", + " 'MarchApril|p_R': array([1.4273317e-12], dtype=float32),\n", + " 'MarchApril|status': array([0], dtype=int32),\n", + " 'MarchApril|urmsd': array([0.03910225], dtype=float32),\n", + " 'MarchApril|urmsd_ci_lower': array([0.03201326], dtype=float32),\n", + " 'MarchApril|urmsd_ci_upper': array([0.05224345], dtype=float32),\n", + " 'gpi': array([0], dtype=int32),\n", + " 'lat': array([33.8833]),\n", + " 'lon': array([102.1333])}}\n", + "Results were stored in /tmp/tmppeisd__e/ascat_ismn_adapted\n" + ] + } + ], + "source": [ + "import pandas as pd\n", + "from pytesmo.validation_framework.temporal_matchers import make_combined_temporal_matcher\n", + "process = Validation(\n", + " datasets, 'ISMN',\n", + " temporal_ref='ASCAT',\n", + " scaling='mean_std',\n", + " scaling_ref='ISMN',\n", + " metrics_calculators={(2, 2): adapted_intercomparison_metrics.calc_metrics},\n", + " temporal_matcher=make_combined_temporal_matcher(pd.Timedelta(6, \"H\")),\n", + " period=period)\n", + "\n", + "save_path = output_folder / 'ascat_ismn_adapted'\n", + "\n", + "if os.path.exists(save_path):\n", + " shutil.rmtree(save_path)\n", + "\n", + "for job in jobs:\n", + " results = process.calc(*job)\n", + " pprint(results)\n", + " netcdf_results_manager(results, save_path)\n", + "\n", + "print(f\"Results were stored in {save_path}\")" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -489,7 +681,7 @@ }, { "cell_type": "code", - "execution_count": 11, + "execution_count": 15, "metadata": {}, "outputs": [], "source": [ @@ -521,27 +713,41 @@ "\n", "setup_code = \"my_validation.py\"\n", "start_validation(setup_code)\n", - "```\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "## Time Series Adapters\n", "\n", - "## Masking datasets\n", + "Data readers extract time series from the candidate and reference data sets used in a validation run. However, often it is desired to change the data in some way after reading and before using them in the validation framework. Potential preprocessing steps include masking/filtering, conversion to anomalies and the combining multiple available variables in a dataset.\n", + "\n", + "## Masking / filtering\n", "\n", - "Masking datasets are datasets that return a pandas DataFrame with boolean values. `True` means that the observation\n", - " should be masked, `False` means it should be kept. All masking datasets are temporally matched in pairs to the\n", - "temporal reference dataset. Only observations for which all masking datasets have a value of `False` are kept for\n", - "further validation.\n", + "Filters are used to select certain observations in a time series after reading to include/exclude from the validation run. Satellite and in situ observations often come with quality flags for their measurements. Filters are for example used to choose - based on these flags - which observations should be used and which ones should be discarded before computing validation metrics.\n", "\n", - "The masking datasets have the same format as the dataset dictionary and can be specified in the Validation class\n", - "with the `masking_datasets` keyword.\n", + "There are 2 ways of masking input datasets.\n", + "\n", + "1) Directly upon loading a time series from a dataset, by removing any unwanted observations\n", + " For this the `SelfMaskingAdapter` and `AdvancedMaskingAdapter` are used.\n", + "2) Indirectly after reading all datasets by using one or multiple independent masking datasets.\n", + " For this the `MaskingAdapter` is used. The masking datasets have the same format as the dataset dictionary and can be specified in the Validation class with the `masking_datasets` keyword.\n", "\n", "### Masking adapter\n", "\n", - "To easily transform an existing dataset into a masking dataset `pytesmo` offers a adapter class that calls the\n", + "To easily transform an existing dataset into a masking dataset `pytesmo` offers an adapter class that calls the\n", "reading method of an existing dataset and creates a masking dataset based on an operator, a given threshold, and (optionally) a column name." ] }, { "cell_type": "code", - "execution_count": 12, + "execution_count": 16, "metadata": {}, "outputs": [ { @@ -554,23 +760,23 @@ "2008-07-01 01:00:00 False\n", "2008-07-01 02:00:00 False\n", "2008-07-01 03:00:00 False\n", - "2008-07-01 04:00:00 False\n" - ] - }, - { - "name": "stderr", - "output_type": "stream", - "text": [ - "/tmp/ipykernel_19692/3230815451.py:3: DeprecationWarning: `MaskingAdapter` is deprecated, use `SelfMaskingAdapter` or `AdvancedMaskingAdapter` instead.\n", - " ds_mask = MaskingAdapter(ismn_reader, '<', 0.2, 'soil_moisture')\n" + "2008-07-01 04:00:00 False\n", + "... ...\n", + "2010-07-31 19:00:00 False\n", + "2010-07-31 20:00:00 False\n", + "2010-07-31 21:00:00 False\n", + "2010-07-31 22:00:00 False\n", + "2010-07-31 23:00:00 False\n", + "\n", + "[15927 rows x 1 columns]\n" ] } ], "source": [ "from pytesmo.validation_framework.adapters import MaskingAdapter\n", "\n", - "ds_mask = MaskingAdapter(ismn_reader, '<', 0.2, 'soil_moisture')\n", - "print(ds_mask.read(ids[0]).head())" + "ds_mask = MaskingAdapter(ismn_reader, '<', 0.2, 'soil_moisture')\n", + "pprint(ds_mask.read(ids[0]))" ] }, { @@ -578,12 +784,12 @@ "metadata": {}, "source": [ "### Self-masking adapter\n", - "`pytesmo` also has a class that masks a dataset \"on-the-fly\", based on one of the columns it contains and an operator and a threshold. In contrast to the masking adapter mentioned above, the output of the self-masking adapter is the masked data, not the the mask. The self-masking adapter wraps a data reader, which must have a `read_ts` or `read` method. Alternatively, to use a method with a name other than `read`/`read_ts`, use the `read_name` keyword which is available for each Adapter - it is still required that the method returns a pandas DataFrame! Calling the reading method will return the masked data - more precisely a DataFrame with only rows where the masking condition is true." + "`pytesmo` also has a class that masks a dataset \"on-the-fly\", based on one of the columns it contains and an operator and a threshold. In contrast to the masking adapter mentioned above, the output of the self-masking adapter is the masked data, not the mask. The self-masking adapter wraps a data reader, which must have a `read_ts` or `read` method. Alternatively, to use a method with a name other than `read`/`read_ts`, use the `read_name` keyword which is available for each Adapter - it is still required that the method returns a pandas DataFrame! Calling the reading method will return the masked data - more precisely a DataFrame with only rows where the masking condition is true." ] }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 17, "metadata": {}, "outputs": [ { @@ -596,16 +802,212 @@ "2008-11-29 10:00:00 0.19 D01,D03 M\n", "2008-11-29 11:00:00 0.19 D01,D03 M\n", "2008-11-30 03:00:00 0.19 D01,D03 M\n", - "2008-11-30 04:00:00 0.19 D01,D03 M\n" + "2008-11-30 04:00:00 0.19 D01,D03 M\n", + "... ... ... ...\n", + "2010-03-16 08:00:00 0.18 D01 M\n", + "2010-03-16 09:00:00 0.18 D01 M\n", + "2010-03-16 10:00:00 0.18 D01 M\n", + "2010-03-16 11:00:00 0.18 D01 M\n", + "2010-03-16 12:00:00 0.19 D01 M\n", + "\n", + "[2956 rows x 3 columns]\n" ] } ], "source": [ "from pytesmo.validation_framework.adapters import SelfMaskingAdapter\n", "\n", - "ds_mask = SelfMaskingAdapter(ismn_reader, '<', 0.2, 'soil_moisture')\n", - "print(ds_mask.read(ids[0]).head())" + "ds_mask = SelfMaskingAdapter(ismn_reader, '<', 0.2, 'soil_moisture', read_name='read')\n", + "pprint(ds_mask.read(ids[0]))" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "While the `SelfMaskingAdapter` works only one one filtering rule, there is also the `AdvancedMaskingAdapter`. This one can take a list of conditions based on which a time series is masked." + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "metadata": { + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " sm snow_prob proc_flag orbit_dir\n", + "2007-02-12 14:47:48.999982 16.0 0 0 b'A'\n", + "2007-02-20 13:42:25.999977 35.0 0 0 b'A'\n", + "2007-03-01 13:56:14.999987 37.0 0 0 b'A'\n", + "2007-03-06 13:52:51.000007 35.0 0 0 b'A'\n", + "2007-03-08 14:51:29.999990 36.0 0 0 b'A'\n", + "... ... ... ... ...\n", + "2014-06-20 13:54:58.000002 80.0 0 0 b'A'\n", + "2014-06-22 14:53:39.000015 77.0 0 0 b'A'\n", + "2014-06-25 13:51:37.000006 83.0 0 0 b'A'\n", + "2014-06-27 14:50:18.000019 81.0 0 0 b'A'\n", + "2014-06-30 13:48:17.999999 83.0 0 0 b'A'\n", + "\n", + "[656 rows x 4 columns]\n" + ] + } + ], + "source": [ + "from pytesmo.validation_framework.adapters import AdvancedMaskingAdapter\n", + "\n", + "\n", + "ds_mask = AdvancedMaskingAdapter(ascat_reader,\n", + " filter_list=[('snow_prob', np.less_equal, 10),\n", + " ('sm', '>', 0),\n", + " ('sm', '<', 100),\n", + " ('proc_flag', '==', 0),\n", + " ('orbit_dir', '==', b'A')],\n", + " read_name='read')\n", + "pprint(ds_mask.read(1814367)[['sm', 'snow_prob', 'proc_flag', 'orbit_dir']])" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "## AnomalyAdapter and AnomalyClimAdapter\n", + "\n", + "These 2 adapters are used to transform absolute values into anomalies on-the-fly after reading. You can select one or multiple columns for which the anomalies are computed. Additional kwargs (like `timespan`) are passed to the underlying anomaly and climatology function. For more details there is a separate tutorial available." + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": { + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " soil_moisture soil_moisture_flag soil_moisture_orig_flag\n", + "date_time \n", + "2008-07-01 00:00:00 0.140861 C03 M\n", + "2008-07-01 01:00:00 0.140861 C03 M\n", + "2008-07-01 02:00:00 0.140861 C03 M\n", + "2008-07-01 03:00:00 0.140861 C03 M\n", + "2008-07-01 04:00:00 0.140861 C03 M\n", + "... ... ... ...\n", + "2010-07-31 19:00:00 -0.151816 U M\n", + "2010-07-31 20:00:00 -0.151816 U M\n", + "2010-07-31 21:00:00 -0.151816 U M\n", + "2010-07-31 22:00:00 -0.151816 U M\n", + "2010-07-31 23:00:00 -0.151816 U M\n", + "\n", + "[15927 rows x 3 columns]\n" + ] + }, + { + "data": { + "text/plain": "" + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "text/plain": "
", + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAnEAAAD2CAYAAAC5p9FxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACOk0lEQVR4nO2dd3wT5R/HP0n3Lt0tdEFZhbI3yJA9lCFLEWTIliGiwE9lCAIOEFSWiCBLQEFQlgzZIKPsXUZpoS0t3TNNk+f3R3qXu+Sy2qRJmuf9evXV5PLc3fMk99x9n+8UEUIIKBQKhUKhUChWhdjcHaBQKBQKhUKhGA4V4igUCoVCoVCsECrEUSgUCoVCoVghVIijUCgUCoVCsUKoEEehUCgUCoVihVAhjkKhUCgUCsUKoUIchUKhUCgUihVChTgKhUKhUCgUK4QKcRQKhUKhUChWCBXiKJQK4Pvvv4dIJEL9+vXN3RWTEBERgZEjRxr1mI8fP4aTkxMuXLhg1ONySUpKwvz583H9+vVyHWfz5s0YOnQoateuDbFYjIiICI1tL126hO7du8PDwwPu7u7o1KkTzp07p9aOEILvv/8ederUgZOTE4KDgzFx4kRkZmaqtV2xYgUGDBiAyMhIiEQidOzY0aD+//vvvxg9ejTq1KkDNzc3VK1aFX379kVsbKxg+6tXr6JLly5wd3eHt7c3BgwYgCdPnvDaPHz4EDNnzkTTpk3h7e0NHx8ftG3bFn/88YfgMVNTUzFy5Ej4+fnB1dUVrVu3xvHjxw0ah6XQvn17TJ8+3dzdoNgAVIijUCqAX375BQBw584dXLx40cy9sQ5mzpyJrl27onXr1iY7R1JSEhYsWFBuIW7Lli24c+cOWrRogRo1amhsd/nyZbRv3x6FhYXYsmULtmzZgqKiInTu3FlNWJ05cyY+/PBD9O3bF/v378fs2bOxfft2dO3aFVKplNd27dq1ePbsGV5//XX4+/sb3P81a9YgPj4e06ZNw8GDB7Fy5UqkpqaiVatW+Pfff3lt79+/j44dO6K4uBi7du3CL7/8gocPH+K1115DWloa2+7IkSM4cOAA3nrrLfz+++/Ytm0batasiUGDBuGLL77gHVMikaBz5844fvw4Vq5ciX379iEwMBA9evTAqVOnDB6PuVm4cCFWr16NBw8emLsrlMoOoVAoJuXy5csEAOnduzcBQMaOHWvuLhmd8PBw8t577xnteHfv3iUAyOHDh412TC4lJSWkqKiI/W02btxYruPJZDL2de/evUl4eLhgu+7du5PAwECSn5/PbsvJySF+fn6kTZs27Lbnz58TOzs7MmXKFN7+27dvJwDITz/9pPH89erVIx06dDCo/y9fvlTblpubSwIDA0nnzp152wcNGkT8/PxIdnY2uy0+Pp44ODiQTz75hN2WlpZG5HK52nF79+5NXF1dSVFREbtt1apVBAA5f/48u00qlZLo6GjSokULg8ZiKdSvX79SznWKZUE1cRSKidmwYQMAYOnSpWjTpg127NiBgoICXpv4+HiIRCJ8++23WL58OSIjI+Hu7o7WrVvjv//+UzvmX3/9hdatW8PV1RUeHh7o2rWrmiZn/vz5EIlEuHnzJgYNGgQvLy/4+PhgxowZKCkpwYMHD9CjRw94eHggIiICX3/9NW//oqIifPTRR2jUqBG7b+vWrbFv3z6t483Ly4O3tzfGjx+v9ll8fDzs7OzwzTffaD3GmjVrEBQUhK5du7LbVq1aBbFYjNTUVHbbsmXLIBKJMHnyZHabXC5HlSpV8NFHH7HnFIlE+Prrr7Fo0SJERkbCyckJJ06cQPPmzQEAo0aNgkgkgkgkwvz587X2TQixWL9b6blz59CxY0e4urqy2zw8PNC+fXucP38eycnJAID//vsPMpkMvXr14u3fp08fAMDu3bvLdH5NBAQEqG1zd3dHdHQ0EhMT2W0lJSXYv38/3nrrLXh6erLbw8PD0alTJ/z555/sNj8/P4hEIrXjtmjRAgUFBcjIyGC3/fnnn6hduzZP62pvb493330Xly5dwosXL7T2/+jRo+jbty+qVasGZ2dnREVFYfz48Xj16hXbZu/evRCJRIIm2jVr1rBzhWH9+vWoVasWnJycEB0dje3bt2PkyJFaTeVchg8fju3btyM3N1ev9hRKWaBCHIViQgoLC/Hbb7+hefPmqF+/PkaPHo3c3Fz8/vvvgu1XrVqFo0ePYsWKFdi2bRvy8/PRq1cvZGdns222b9+Ovn37wtPTE7/99hs2bNiAzMxMdOzYEWfPnlU75uDBg9GwYUPs3r0bY8eOxXfffYcPP/wQ/fr1Q+/evfHnn3/i9ddfx6xZs7Bnzx52P4lEgoyMDMycORN79+7Fb7/9hnbt2mHAgAHYvHmzxjG7u7tj9OjR2LZtG6/fALB69Wo4Ojpi9OjRWr+3AwcOoH379jzhpEuXLiCE8B7Cx44dg4uLC44ePcpuu3LlCrKystClSxfeMb///nv8+++/+Pbbb3Ho0CGEhIRg48aNAIDPPvsMFy5cwIULF/D+++9r7Vt5KC4uhpOTk9p2ZtutW7fYdtztDA4ODmrChqnIzs7G1atXUa9ePXbb48ePUVhYiAYNGqi1b9CgAR49eoSioiKtxz1x4gT8/f15guPt27c1HhNQuCFo4/Hjx2jdujXWrFmDI0eOYO7cubh48SLatWvHmp779OmDgIAA9jfnsmnTJjRp0oQ9308//YRx48ahQYMG2LNnDz777DMsWLAAJ0+e1NoPLh07dkR+fr5B+1AoBmNuVSCFUpnZvHkzAUDWrl1LCFGYqNzd3clrr73Ga/f06VMCgMTExJCSkhJ2+6VLlwgA8ttvvxFCFGazkJAQEhMTwzOh5ebmkoCAAJ5Jbt68eQQAWbZsGe9cjRo1IgDInj172G1SqZT4+/uTAQMGaBxLSUkJkUqlZMyYMaRx48a8z1TNqY8fPyZisZh899137LbCwkLi6+tLRo0apfEchChMewDI0qVL1T6rVq0aGT16NCGEEIlEQtzc3MisWbMIAPLs2TNCCCFffvklcXBwIHl5eYQQ5Xdbo0YNUlxczDuescypXLSZUxs1akRq1arF++2kUimpXr06AUC2b99OCCHk+vXrBABZuHAhb//jx48TAMTR0VHj+ctiThVi2LBhxN7enly5coXddu7cOd71yGXx4sUEAElKStJ4zPXr1xMAZOXKlbztDg4OZPz48Wrtz58/z/te9EEulxOpVEqePXtGAJB9+/axn82YMYO4uLiQrKwsdhtjuv/hhx8IIYo5FhQURFq2bMk77rNnz4iDg4PG31aV4uJiIhKJyKxZs/TuO4ViKFQTR6GYkA0bNsDFxQVDhw4FoNBSDRo0CGfOnEFcXJxa+969e8POzo59z2gGnj17BgB48OABkpKSMHz4cJ6Wyt3dHW+99Rb+++8/NVMtY4JjqFu3LkQiEXr27Mlus7e3R1RUFHseht9//x1t27aFu7s77O3t4eDggA0bNuDevXtax129enX06dMHq1evBiEEgEKDmJ6ejg8++EDrvklJSQCETXydO3fGsWPHAADnz59HQUEBZsyYAT8/P1Ybd+zYMbRu3Rpubm68fd988004ODhoPbepmTJlCh4+fIgPPvgAL168QGJiIiZMmMB+78xv2rBhQ7Rv3x7ffPMNfv/9d2RlZeH8+fOYMGEC7OzsymQ+JYSgpKSE96eJzz//HNu2bcN3332Hpk2bqn0uZCbV9dmhQ4cwefJkDBw4EFOmTDHKMRlSU1MxYcIEhIaGstdpeHg4APCu1dGjR6OwsBA7d+5kt23cuBFOTk545513ACjmWEpKCgYPHsw7R1hYGNq2bau1H1wcHBzg7e2t0xRMoZQHKsRRKCbi0aNHOH36NHr37g1CCLKyspCVlYWBAwcCUEascvH19eW9Z8xphYWFAID09HQAQHBwsNq+ISEhkMvlaikofHx8eO8dHR3h6uoKZ2dnte1cU9iePXswePBgVK1aFVu3bsWFCxdw+fJljB49WqfJDACmTZuGuLg4VrhatWoVWrdujSZNmmjdjxmrav8AhUk1ISEBcXFxOHbsGBo3boyAgAC8/vrrOHbsGAoLC3H+/Hk1Uyog/J1VNKNHj8bSpUuxZcsWVKtWDWFhYbh79y5mzpwJAKhatSrblhGgBw8ejCpVqqBTp04YMGAAGjVqxGunL6dOnYKDgwPvLz4+Xq3dggULsGjRInz55ZdqAjdzfTLXIZeMjAyIRCJ4e3urffbPP/9gwIAB6Nq1K7Zt26YmlPn6+mo8JqB+DXORy+Xo1q0b9uzZg08++QTHjx/HpUuXWF9S5noCgHr16qF58+asSVUmk2Hr1q3o27cvew6mH4GBgWrnEtqmDWdnZ975KRRjY2/uDlAolZVffvkFhBD88ccfgrmxfv31VyxatIinedMF8xBlHOC5JCUlQSwWo0qVKmXvNIetW7ciMjISO3fu5D10JRKJXvu//vrrqF+/Pn788Ue4u7vj6tWr2Lp1q879/Pz8AIDn+M7QuXNnAApt29GjR9nAh86dO+Ozzz7D6dOnIZFIBIU4XdqcimLWrFmYPn064uLi4OHhgfDwcIwfPx5ubm48rVdAQAAOHjyI1NRUpKSkIDw8HC4uLli9ejW7EDCEpk2b4vLly7xtISEhvPcLFizA/PnzMX/+fPzvf/9TO0aNGjXg4uLC+u5xuXXrFqKiotSE73/++Qf9+vVDhw4dsHv3bjg6OqrtGxMTo/GYALTmV7x9+zZu3LiBTZs24b333mO3P3r0SLD9qFGjMGnSJNy7dw9PnjxBcnIyRo0axX7OzLGXL1+q7ZuSkqKxH0JkZmay1zOFYgqoJo5CMQEymQy//voratSogRMnTqj9ffTRR0hOTsahQ4cMOm7t2rVRtWpVbN++nTVTAkB+fj52797NRqwaA5FIBEdHR57wk5KSojM6lcvUqVNx4MABzJkzB4GBgRg0aJDOfRhh5fHjx2qfBQcHIzo6Grt370ZsbCwrxHXt2hVpaWlYvnw5PD092ahTXahqOisKJycn1K9fH+Hh4UhISMDOnTsxduxYuLi4qLUNCAhAgwYN4OXlhbVr1yI/P1+nSVoIDw8PNGvWjPfHFagWLlyI+fPn47PPPsO8efMEj2Fvb4833ngDe/bs4UVdJiQk4MSJExgwYACv/ZEjR9CvXz+0a9cOe/fuFQzqAID+/fvj/v37vByKJSUl2Lp1K1q2bKkmbHJhrk/VY69bt06w/dtvvw1nZ2ds2rQJmzZtQtWqVdGtWzf289q1ayMoKAi7du3i7ZeQkIDz589r7IcqSUlJKCoqQnR0tN77UCiGQjVxFIoJOHToEJKSkvDVV18JZs9nNFQbNmxQ81nThlgsxtdff41hw4ahT58+GD9+PCQSCb755htkZWVh6dKlRhtDnz59sGfPHkyaNAkDBw5EYmIiFi5ciODgYEF/PiHeffddzJkzB6dPn8Znn30mqIVRxdHRUWNqFUChdfvhhx/g4uLC+ihFRkYiMjISR44cwZtvvgl7e/1ubYxmadu2bahbty7c3d0REhKCkJAQbN68GaNHj8Yvv/yCESNGaD3O3bt3cffuXQAKQbegoIDVvkZHR7MP8tu3b2P37t1o1qwZnJyccOPGDSxduhQ1a9bEwoULecdcv34928esrCwcOnQIGzZswOLFi9VM0leuXGFNozk5OawGGACaN2/O+odpYtmyZZg7dy569OiB3r17q333rVq1Yl8vWLAAzZs3R58+fTB79mwUFRVh7ty58PPzY9O6AMDZs2fRr18/BAUF4X//+59aQuXo6Gg2Tcno0aOxatUqDBo0CEuXLkVAQACbLJfxgdREnTp1UKNGDcyePRuEEPj4+ODvv//mRSxz8fb2Rv/+/bFp0yZkZWVh5syZPB9DsViMBQsWYPz48Rg4cCBGjx6NrKwsLFiwAMHBwWr+iFFRUQDUNX/Md9ipUyet/adQyoUZgyoolEpLv379iKOjI0lNTdXYZujQocTe3p6kpKSwEZTffPONWjsAZN68ebxte/fuJS1btiTOzs7Ezc2NdO7cmZw7d47XholOTUtL421/7733iJubm9p5OnToQOrVq8fbtnTpUhIREUGcnJxI3bp1yfr169njctGW7HfkyJHE3t6ePH/+XNNXocaGDRuInZ2dYKTjvn37CADStWtX3vaxY8cSAOT777/nbdf23RJCyG+//Ubq1KlDHBwceN/1xo0b9Y5cZb4ToT/ub/fgwQPSvn174uPjQxwdHUlUVBT57LPP2EhaLuvWrSN169Ylrq6ubETz3r17Bc//3nvvaTy/Pv3v0KGDxv2FHhNXrlwhnTt3Jq6ursTT05P069ePPHr0SO/vBAA5ceIEr31KSgoZMWIE8fHxIc7OzqRVq1bk6NGjOvtOiCLCtGvXrsTDw4NUqVKFDBo0iCQkJAjOHUIIOXLkCNuPhw8fCh7zp59+IlFRUcTR0ZHUqlWL/PLLL6Rv376CkdlCEavDhw8nMTExevWfQikrIkI4NhkKhUIxIsXFxYiIiEC7du3UzFPaKCoqQlhYGD766CPMmjXLhD2kUPQjKysLtWrVQr9+/fDTTz9pbZuTk4OQkBB89913GDt2bAX1kGKLUJ84CoVidNLS0nD27FlMnDgRL1++xOzZsw3a39nZGQsWLMDy5cuRn59vol5SKMKkpKRgypQp2LNnD06dOoXNmzejU6dOyM3NxbRp03Tu/9133yEsLIwXMEGhmALqE0ehUIzOgQMHMGrUKAQHB2P16tU604oIMW7cOGRlZeHJkyeIiYkxQS8pFGGcnJwQHx+PSZMmISMjA66urmjVqhXWrl3Lq2ChCU9PT2zatElv30wKpaxQcyqFQqFQKBSKFULNqRQKhUKhUChWiNUJcatXr0ZkZCScnZ3RtGlTnDlzRmPbPXv2oGvXrvD394enpydat26Nf/75pwJ7S6FQKBQKhWIarMqcunPnTgwfPhyrV69G27ZtsW7dOvz888+4e/cuwsLC1NpPnz4dISEh6NSpE7y9vbFx40Z8++23uHjxIho3bqzXOeVyOZKSkuDh4WExGd8pFAqFQqFUXgghyM3NRUhIiNZayVYlxLVs2RJNmjTBmjVr2G1169ZFv379sGTJEr2OUa9ePQwZMgRz587Vq/3z588RGhpapv5SKBQKhUKhlJXExERUq1ZN4+dWEzpTXFyM2NhYtVQF3bp107sUilwuR25urtZiyhKJhFcbkpFxExMT2eziFAqFQqFQKKYiJycHoaGh8PDw0NrOaoS4V69eQSaTITAwkLc9MDBQ76LEy5YtQ35+PgYPHqyxzZIlS7BgwQK17Z6enlSIo1AoFAqFUmHocuOyusAG1QERQvTyVfvtt98wf/587Ny5EwEBARrbzZkzB9nZ2exfYmJiuftMoVAoFAqFYmysRhPn5+cHOzs7Na1bamqqmnZOlZ07d2LMmDH4/fff0aVLF61tnZyc4OTkVO7+UigUCoVCoZgSq9HEOTo6omnTpjh69Chv+9GjR9GmTRuN+/32228YOXIktm/fjt69e5u6mxQKpYJIzCjAXzeSIJdbTWwWhUKhGBWr0cQBwIwZMzB8+HA0a9YMrVu3xk8//YSEhARMmDABgMIU+uLFC2zevBmAQoAbMWIEVq5ciVatWrFaPBcXF3h5eZltHBQKpfy89vUJAIBEKsOgZjSCnGIZyGQySKVSc3eDooOHKbl4kVWITnU0u1eZEgcHB9jZ2ZX7OFYlxA0ZMgTp6en44osvkJycjPr16+PgwYMIDw8HACQnJyMhIYFtv27dOpSUlGDy5MmYPHkyu/29997Dpk2bKrr7ZUIuJ8iVlMDLxcHcXaFQLJKLTzOoEEcxO4QQpKSkICsry9xdoejBy8xC2AN4EJcLR3vzGCW9vb0RFBRUrhy0VpUnzhzk5OTAy8sL2dnZZolOjZh9AACw9t0m6FE/uMLPT7Ft8iQluJ6QhVbVfWBvZ1neF8zcGNi0Gr4d1NDMvaHYOsnJycjKykJAQABcXV1pcngL50FKDgAg0NMZ3q6OFXpuQggKCgqQmpoKb29vBAerP9v1lT2sShNny0zYehXxS6lPH6ViWfj3Xey8koiF/epjeKtwc3eHQrFIZDIZK8D5+vqauzsUPRDZFwEAHJ2c4OysCGYsKC4BCODqZHrRyMXFBYAiODMgIKDMplXLWlpTKBSLYucVRYqdtScfm7knmqG2BIq5YXzgXF1dzdwTSlmRE4JHqXl4lJYHWQUFSzHXS3l8KKkQR6FQdCK3YEmJwHL7RrEtqAnVeuHe4ypKiDPG9UKFOAqFohNLFuKoDEehUMoKc/uwVvGbCnEUCkUnlizDvcwtQlJWobm7QaFQrBELvrfpAxXiKBSKTiz5PnfuUTraLP0XhcUyc3eFQrEZ4uPjIRKJcP36dQDAyZMnIRKJzJpiZdOmTfD29i7j3lxdnCXf8fhQIY5CoejEGjIRpeVKzN0FCsVmCA0NZfO1WgpDhgzBw4cP9Wq7b9d2tKtn/RH3NMWIldAsvIq5u0CxYayhspVF++1RKJUMOzs7BAUFmbsbPFxcXNjUHeXB0DuJVCqFg4N5EvJTTZyV4EkrNlDMSEZ+sbm7oBMqwlEsCUIICopLKvzPUK35H3/8gZiYGLi4uMDX1xddunRBfn4+5HI5vvjiC1SrVg1OTk5o1KgRDh8+zO6nak41BMbsuX//ftSuXRuurq4YOHAg8vPz8euvvyIiIgJVqlTBlClTIJMp3SQyMzMxYsQIVKlSBa6urujZsyfi4uLUjstw48YNdOrUCR4eHvD09ETTpk1x5coVnDx5EnM/mozcnBz4eThDJBJhwYL5AICGoVXw1759vP56e3uzVZ6Yce/atQsdO3aEs7Mztm7dCgDYuHEj6tatC2dnZ9SpUwerV682+LsxFKqJsxKkMrm5u0CxcZ6+ykekn5u5u6ERqomjWBKFUhmi5/5T4ee9+0V3uDrq92hPTk7G22+/ja+//hr9+/dHbm4uzpw5A0IIVq5ciWXLlmHdunVo3LgxfvnlF7z55pu4c+cOatasWe5+FhQU4Pvvv8eOHTuQm5uLAQMGYMCAAfD29sbBgwfx5MkTvPXWW2jXrh2GDBkCABg5ciTi4uLw119/wdPTE7NmzUKvXr1w9+5dQU3YsGHD0LhxY/zw4yrkSwkuXo6FSGyHNm3a4JP5S7B62WLsO3kZ0SFecHF1w7OcUoFRj1vJrFmzsGzZMmzcuBFOTk5Yv3495s2bhx9//BGNGzfGtWvXMHbsWLi5ueG9994r9/elCSrEWTDFJUrB7UZilvk6QqEAyJeUmLsLWqEyHIViGMnJySgpKcGAAQPYGuQxMTEAgG+//RazZs3C0KFDAQBfffUVTpw4gRUrVmDVqlXlPrdUKsWaNWtQo0YNAMDAgQOxZcsWvHz5Eu7u7oiOjkanTp1w4sQJDBkyhBXezp07hzZt2gAAtm3bhtDQUOzduxeDBg1SO0dCQgJGT5wKuVcIXAB07FkVfu5OcHR0hLuHJ0QiEfwCAhEU5K3IDZeTrXf/p0+fjgEDBrDvFy5ciGXLlrHbIiMjcffuXaxbt44KcbZKiVwpxMVU8zJjTygUy8cagi8otoOLgx3uftHdLOfVl4YNG6Jz586IiYlB9+7d0a1bNwwcOBB2dnZISkpC27Ztee3btm2LGzduGKWfrq6urAAHAIGBgYiIiIC7uztvW2pqKgDg3r17sLe3R8uWLdnPfX19Ubt2bdy7d0/wHDNmzMCs6ZPRdMd2tGzXAd1694NPdG2j9L9Zs2bs67S0NCQmJmLMmDEYO3Ysu72kpAReXqZ9dlMhzswkZxfiblIO6oV4IcjLmfcZ15lcIqXmVErFkl3ALwVTYuHRDZbdO4qtIRKJ9DZrmgs7OzscPXoU58+fx5EjR/DDDz/g008/xdGjRwGoVxQghBitKoWq+VMkEgluk5cqMzQt0rT1af78+Wjyeh+cOX4EZ08cw5rlS7Fq/a8YO2Koxn6JRCI11wyhslhubkrXEqaP69ev5wmZAMpcE1VfaGCDmTn/KB1jfr2CPj+chVzlIcm9aK88y0SRlObBolQcu0rrpjI8SMkx6fmepOXh87238aKMiXupTxyFYjgikQht27bFggULcO3aNTg6OuL48eMICQnB2bNneW3Pnz+PunXrmqWf0dHRKCkpwcWLF9lt6enpePjwodY+RVSPwvCxk7Bu+x507tEHu37bAkAhRMo4vubM87aKrx9SUpLZ7XFxcSgoKNDat8DAQFStWhVPnjxBVFQU7y8yMrJM49UXy14m2ACujgop/VWeBJISOVwclVK7quLjw53XsebdpmU6jzFXUBTbQPVymbX7FhqFVkHtIA+TnG/oT/8hNVeCK88ycWjaa1rbxr3MVdtGZTgKxTAuXryI48ePo1u3bggICMDFixeRlpaGunXr4uOPP8a8efNQo0YNNGrUCBs3bsT169exbds2s/S1Zs2a6Nu3L8aOHYt169bBw8MDs2fPRtWqVdG3b1+19oWFhfj444/RuEMPVA0Nw8vkJNy5cQ293+wHmZwgJDQMBfl5uHj2FEK6tEVSngywd0KLNq9h7erVeK1tG8jlcsyaNUuv9CHz58/H1KlT4enpiZ49e0IikeDKlSvIzMzEjBkzTPGVAKBCnNnpEh3Ivi6WyeECpRCnqj4+dDulTOeImH0AADCxYw3M6lGnTMegUACg36pzuLewh0mOnVqarPdesm6N33WBQJ+KKlpNoVQWPD09cfr0aaxYsQI5OTkIDw/HsmXL0LNnT3Tv3h05OTn46KOPkJqaiujoaPz1119GiUwtKxs3bsS0adPQp08fFBcXo3379jh48KCgkGVnZ4f09HR8Nn0C0l+lwbuKLzr37IMZsz9DYkYBGjVriUHvjsInk0ZjXGYGJnw4CxNnzMZHny/CV/+bhvbt2yMkJAQrV65EbGyszr69//77cHV1xTfffINPPvkEbm5uiImJwfTp003wTSgREeoNrJWcnBx4eXkhOzsbnp6eRj8+IQSRcw4CAK581gV+7k7sZ5n5xWi88CivffzS3gafgxHiyro/xTZZf/oJvjyo7jBsqmvIkOt0+8UE/O/PW7xtm0e3QPta/ibpG4WijaKiIjx9+hSRkZFwdnbWvQOlwrj5PIv3voqrIzIL+HkvG1Tz5rWrFegBZwMCRMqKtutGX9mD+sSZGZFIBEc7xc+gmguO+vhQzIm3q+UmmH6Umqe2jWriKBSKLvS5S1jTnYQKcRaAg53C+UhaohLYINCWKk4pFYWTkVaiqblFOHInxahC1svcIrVtNPCHQrEMevbsCXd3d8G/xYsXm7t7urGixyz1ibMApKUPt1f5EoT5urLbhTRxzzMLEerjqradQjE2zvbCazy5nEAs1j9Ipvf3Z5GWK8HCvvUwvHWEUfrGaK+5JGRojyCjUCgVw88//4zCQuEocx8fnwrujQoCAlphseoC0HqkOCrEWQBMZYYBq8/zfYEErqNiWn6LYkLS8yTYdeU53mpSVaMm7lpiFpqGV9H7mGmlAQuH76RoFeLEIvWIbCGKS+RIyVbXxFmy+ZdCsSWqVq1q7i4A0GS5Ut8Wl5qro4XlQs2pFozQA43WUKWYkqk7ruGrw/fx3sbLrJlflfFbrpTp2NwyckL4uDmyr1UTDXMZtekSLjxJV9teRBNiU8wMdXexLJ6lq2vn9fKJq6Cf0RjXCxXiLBgicLmp+s1RKMbk3COFcHQvOUfj3e5VXrHwBzq4HJ+J7ELNwlnL6r7s6/R8ic4+qpJn4bVdKZUXJsWFrqSwlIolp0jz/UYbFfWUZa4XffLQaYKaUy0YRhPnaC9GkKczEjIKDDanWsLKkBCCQ7dTEFPVi/rzWRHcK2dSxxpYffKx1vaJGQUI9nKGvYC/GsP4LVewY1xrwc/sOX52Upnh1+2uK4mY3CnK4P0olPJiZ2cHb29vts6nq6srTa5uAZAS9QVnSbEckMm0PhslRYWwJ6ZzzyCEoKCgAKmpqfD29i5XaS4qxFkAYT6ugk7ZzEUmFnEiWA0W4srfv/Ly981kTP3tGgCap87S0eSX5uGs/YZ27O5LvL/5CjrW9semUS00truRmA25nOBSfAZqB3qgCseEyk3gWxa3gUBPmp+LYj6CgoIAgBXkKMYjX1KCzAIpPF3s4anjXsQlNVM9uCLXwQ5FJTKtz0Z5jmOF5Inz9vZmr5uyYpAQl56ejps3b6Jhw4bw8fHBq1evsGHDBkgkEgwaNKhCaqqtXr0a33zzDZKTk1GvXj2sWLECr70mXKInOTkZH330EWJjYxEXF4epU6dixYoVJu+joSwdEIN3fr6IqAB33nbmIhNBBIdS7caz9Hy04piddKEp19yJ+6k49TAN/+tVF44aohCNxX8C/ksUy8ReLGa1vdxLp3UNX9Twd8PjtHwAgKREBid75U3ul3NPAQAnH6RpPX6hVIYjd1MwYetVtKnhi+1jW7GfVfdzY31YyhLAU0L9RSlmRCQSITg4GAEBAYIF0yllp/Oyk+zr4x911Hu/9/ecVNvWNsoPV59lolBLSqKFfeujTaSfAT00HAcHh3Jp4Bj0FuIuXbqEbt26IScnB97e3jh69CgGDRoEe3t7EEKwdOlSnD17Fk2aNCl3pzSxc+dOTJ8+HatXr0bbtm2xbt069OzZE3fv3kVYWJhae4lEAn9/f3z66af47rvvTNav8sLUS1VNYCrnaOJyixT+PrN238KQ5upj1YSqVkUmJ7ATizBq02UAQISvK0a2NW2BXksw6VL0QywGUHpv4/pkNgr1xvLBjdB31TkAwNb/EjCmnfK6ERtgOvr7pqK49PnHfOGeG5gg1REEIQRN9kuxBOzs7IzycKYoeZGrFLgMqYjRIioIf157wduWWiBHaoEcOUWahbiUfLnVVN7QWwXz6aefYtCgQcjOzsb//vc/9OvXD507d8bDhw8RFxeHd955BwsXLjRlX7F8+XKMGTMG77//PurWrYsVK1YgNDQUa9asEWwfERGBlStXYsSIEfDy8tLrHBKJBDk5Obw/U+PupJSlzz16xb5mNXEiERqGKvt/63m23sdW1cSpJkR9kSWcy8eYUBnOeuAKY8zvVjdYUfIlxNuF/Wzh/ruQc4QmbTKcp7Py+h7dNpKX440r4HMjTvXVxC3sVx8Dm1YDUDY/OgqFUnkRirBPSC/Q6rcLADIremjpLcTFxsZixowZ8PDwwLRp05CUlISxY8eyn0+ePBmXL182SScBoLi4GLGxsejWrRtve7du3XD+/HmjnWfJkiXw8vJi/0JDQ412bE3U8FeaUePT89nXjAAmEgHd6ynt5vdT9BcsVa9FVRVyRVig7iSZXhCmGAc7rhBX+p/ZwmiMGfKKS9i0IVzhL0ElrD+A46vm7mwP7m115+VEwX4s2q9es1WIhtW80L+xIidVidx6zanXEjLReslx7L+ZZO6uUCiVBgcBYS3Q0xl2upKVV0Yhrri4GC4uipW4g4MDXF1d4eentBn7+voiPd10vk+vXr2CTCZDYGAgb3tgYCBSUlKMdp45c+YgOzub/UtMFH7IGBOxWIS+jUIA8DNHcx+i3EinEgPMRqqaONUAClPXZy0oLsGtF/prDinmhatRI5xFBKBewWHDmaeo9dkhzNh1HaceKn3hlhziC2Bcjd3PZ57wIpQP3haeu/qmC7EXi9moVm3pSyydiVuvIjm7CB9sVwQAZRUUo9fKM/jptPaIYIrtsPJYHNp99S+bPJuiGyEhrlgm50XCC2FNWn29hbjQ0FA8efKEfb9jxw4EBwez75OTk3lCnalQDdsmhBg1lNvJyQmenp68v4rApTQSZtcVpdDIRqeKRXB3UmpBDIncUxXSVB90phbicgpp7i5rQqicFjO9VE0QK4/HAQD2XOX7nOSrlLDhLjoKimW8VTDXvM81fWiqg6rqX2lvJ2L7/DJHgqyCsuWwMzeqc/qn009wNzkHiw/eN1OPKJbGd8ce4nlmIVadeGTurlg1sc8ydfrwWlNSfb2FuKFDh/JCp3v37s1q5gDgr7/+QosWmlMLlBc/Pz/Y2dmpad1SU1PVtHPWyKs8xerq4cs89kGkjE4FOtQKYNvqynzPRVVpp+owbixn8NUnH+HNH88iVyW5Ik2VZF2IBcypXCZ3qqHzGM84LgGA+jVWwBHy7icrTe3cZun5xYIBMarHsheLePOh2aJjVhmlqroQlZQhsKMsSEpkPE0pxfKhATz6oymoTpcvOLPwzCooxiOVklyWht5C3Lx58zB06FCNn3/66afYvn27UTolhKOjI5o2bYqjR4/yth89ehRt2rQx2Xkriqocp3FGXc7MVbFIBDuxqEwO3Ol5fNW76r5PX/EfuGXl68MPcPN5NjZfeMbbTmU464K3QuWkuGGY1FF3Mt1n6QU8v01VbS9X6HqjYQj7WvXhFDnnIBYf5JtmVR2O7cVinmavRE6wRkdSYktEVQGqy2XHGOQWSREz/wgGrbtg+pNRjIZQJZ/KzNm4V7obaaCs3xSjiWvx5XF0WX4aD19ariBXrgRh586dg0SiEBJcXV3h5ORklE5pYsaMGfj555/xyy+/4N69e/jwww+RkJCACRMmAFD4s40YMYK3z/Xr13H9+nXk5eUhLS0N169fx927d03az7LATSzIXHjMQ415rjL2fUNUvf1X84M+VPe99DTDwJ5qR5MZzBKRlMjwICWXpkDhIObJcPzrDwDcnPTLStRjxRn2tapwxjXpM0mENWmDfjr9hPde9VgikXp6ky3/8RcS1oDqGAzxey0rZ+JeobhEjthnmSY/F8V42NLtKqdIinc3XCzz/sx31TaKn1vVSUdu1PvJCqGNiZIvjyBpaspVsaFnz564fv06qlevbqz+aGXIkCFIT0/HF198geTkZNSvXx8HDx5EeHg4AIVfXkJCAm+fxo0bs69jY2Oxfft2hIeHIz4+vkL6rDece/iDlFzUCvTAimMPAShrVTI+Qzc4me11oeoDp5q6Qcjxszyo3WAsWBU3fMMlXHqage/fbow3ORohW4ar1eKa88uDqiYuT6K8JteeeozpXWpixq7reh1LyJQU6uPCe29NCwkGVc2bu57CcvnOacGTk6IRG5LhUFSsPpeZXKf6wCxE64V48Wou63JXcHXiR+JbcuR7uZ7g5tBgTJo0CfHx8ZBIJIiNjUX79u3ZzzZt2oSTJ0/y2hNC1P4sToAD0KGmP/uayRV35O5LXhvGb648FRBSsot47yUlxn3gqT6wRRYsxTFayO0XrU9zYyqE8sSV17FRVfDKl/CvuW//eYCDt/SLMD+qMicIAYK9+EJcTpH1BdOo+sTZi01bReX4vZe8dEYU68GWNHFC+dyepefj4ctcbL4Qr9P/lfmuXB3t8HYL/dOFqVoGKkIzXlZMe6eg6E2bKD+8VlMR3asa3cdwPyVX6+eqCGktlh99yHtv7GtT9Xianv/mNmFyJ7++qzpbQEh2KOu3k1Ma5KJ6HZ59xDdN/Hz2Ke894/vJcIwjuP3vz1vs69dq+rFaOG71CEM4/+gVWi4+xjuHKciXlGh1kFb93o2sIOcR+ywTY369gqWHlJGv2QXWm57F9rBcgcLYCD0nCqUydPvuNObuu4PtlxIE9uLsX/pfBBGWDGiAte/qV1EqS8WCVWLBKUfKdatYt25dpYgMtRR6xShSthRqENKyODfaawmZalGAquQXV7xGQqaidhYSAK4lZKLF4uPYd/2FwKf6kVVQzCuYbihcX8BXudaZlsIU6IpONYSfzyiEM0MXCp90r817P33ndfY1tzTXljEtWQ2WSxmLVQ//5RJe5kjw/uYrZdpfX7ouP4Uuy0/jSrywD6qqadPOhJq4BynqwqS2OpIUy4K5xVpTGoyyInTv6P39Wfb19YQsrfsrqx4p/uvrPqRaA9qSI97Ldad455134ObmZqy+2DzMg+jYPWGtADfZav/V59GHczELYQ5ll6raWdVMlC8pwcStV5GWK8G0HdfLfJ6O355Ev1XncCZOe8F1TWRyBOJGod5l7kdlg29OVQ9s0MWxGR3Y14WliwhGE9c0vIpex/BydUD80t7K4+ghYDg7CN/Kbr/Ixqd/3kJytnBKgYpK15BU6sZw4JaibuzN51m8EnvqQpzytbG11r+ce6q7EcViISDYd/0Fan56CAdK6xBXVnRF4urOc1p6Dyt9pyrEvdWkmkZLTEa+cnEvrUzm1KKiInzzzTfo1asXmjVrhiZNmvD+KGXHsVRI46Yb4eLt6sh7nysp0boa03Tz//KA6aJzdT0U5/91xygh8oxWsqxmMK5Jz9ZC9rXBj05VYIg5lStM3UhUVOpg0oIMblZNcB9V7FQEGn0ELUaLrcqMXdex7WICz3RoTpixvPnjOQz7+SLro6r6MOJOa00phbILpOi87CS+U3GR0MWj1Dy1bbag1aksyAnYBfDk7VfN2xkTo0tGO1GqMfsj9jlOPEhV+1yXJi7c1xX3F/YQPDa3upHQnLEUDA6BGj16NI4ePYqBAweiRYsWRq2WYOtE+im0mpoiZ34a0RTtvjrB21YolWlUEWuaAOvPmG4lruuBe/BWspowWh7KKn5xv7JcK3SENxVCgQ2qc3xU2whsPBfPvq/i6sBqNt0clbeUS/EZyCooZp2EXRz5t5vq/m54kqbuEqDtnlIvxBN3knLwTsswlWO5C7Z/+FJx8+WWBTMnJXLCW1y9yCpAkJczvFwceO2evlI+NPIkJfCxV58zm87H43FaPlYej8OHXWuVu18U68CWAht0jdXBTownaXmY+fsNAOBp8Ln7M/cUR3v1VD7c56eHkz38PZ3wJC2f59bkrTI/LQmDhbgDBw7g4MGDaNu2rSn6Y9MwuWskGsxH1aq4orqfG55wEvQWSGTwdBa+wLir+0961MbXhx+UuW+5RVK4O9nrFNpVhThVbWB+sQz+HurHIITgwctcRPm7C0YkaSI9Tz9/tpTsIrzMKULDUtNpYbFSUD6koX6nLeLM8y3jmyIYhjQP5Qlxm0e3xBs/Kkz73q78azEjv5gVEJqEefM+axJWRVCIY7SBvWKCcPBWClvgHlCuqLvWNcwXN19SghP3U9E4zJu3iAj2ckaySsS2KSmRyXkCE/MywIOfYzOIE3FbUFwCHzd1IU418XF5+0WxDnZffW7uLlQYusylxSUyvMzRXEtW1cqiqvBQve7bRvnhdKmLzm1OzW9LXuQYbE6tWrUqPDw8TNEXm4d5gOZKSpAvKWEFDqE2DK9UKjIQQhAx+wAiZh/g3eQntNddLkkTt19kI2b+EXxUutrRhj4Xu5CQtuHsU/RYcQYf/3HToL456kjayNBqyXH0XXUO8/+6AwC49SLLoPPYCkL5yVTl9kAPZ977mGpeWNSvPjaNaq4m5HOFehcHO14AQoi3C8J9XQXOpzhGy0hFgk5u3jfmnmtoRLFURjBq02U181NF5GPjUiSV88qOMVpK7rw+fDuZV9VCkzlV1excHiz5IUWxXXRdlTlFJSiUKi0p03Zc4++vYk6tVoV/v2HcCKqULj671Qtk5+eXnGoxljw/DBbili1bhlmzZuHZM5pby9j4uitX21eeZaJaFcVq/PM+0ex2JxUHbtU8b88zlQ7c3AhMocLmDLpW4WtOKcoYqRY6F0I1v47QpS9k/v2xtKjzn9cMi1g11Ol70/l4fLb3llq5MXOnPLEU7DlF6DXdt1wclQLHwn71AQDvtgpHx9oBam27fneafW0nFvGCFOxEIjxLV/qdrB7WBCdmdlSep1SwOXQ7BQVskIScPZYqi/vHsK+vJQhXIeAm/ASAOI6viyFJtMtKoVSG45zApZQchRaQK8RN2HoVa08pS4dpyuem26lbfyw5hQLFdtGnru/HvysX/vuuJ/EWQNwUIwDUNNrMAumPiW0w/41ojb61qlkXLAmDhbhmzZqhqKgI1atXh4eHB3x8fHh/lLLjZG+HGv4Kv7jCYhl7BdpzHlj5Er7/FjflAsDXmjC+XrqUFlGfHtKa9NfZXv/0DaorFqHnjIOdeofKmj3+oh5lw/JUvrOt/yWoaTcseKFVoXAFbCaxrmrCZm7JmuYR+kWcAoCrik+cqiwfU9WL9QsFgMYc8yvjWMxo9oSEOK6f3LXS1AOqZkpVuJq4xMwCLS2NQ+yzTN6cZeaHm6PmOfYqV9hc5OFsPC2iJWekp9gu+qxT0vP5LjXca1lVE6cKo4mr4e+OkW0j1SxdynaW+4Aw+C7w9ttv48WLF1i8eDECAwNpYIORCfZyweO0fOQUSQVrV/asH4yHL+PY96pltbgPYSbyTZ/f6Pcrz/Fuq3DBzzSlbxDiXnKOzjZCmriy5ttV9Wc6E5eG5KwiDG6uzM697IhuX0BDSrlUZrgLBtb3UuVrEYlEmNG1FlJyilA7UN214s6C7qg37x+17aqmb7FYhNqBHnhQWlzaXkW4rxnogWpVXPA8s5AVfLQJcQAwqGk1/B77HEWlixIXLcKRYizK1xWhjfJ3d+K5OTBaA21+oJpcBuy1XK8Z+cX45I8bGNQsFN3rBenslyWbiyi2S1kyBxSXyMG4vQrt7+Jgx1oE9J3zllzKz2Ah7vz587hw4QIaNmxoiv7YPMzD6fLTDMHalaoPpZXH4ngqYO7DbeVxhbCnT4oGVWFQqE/64K6iHRCaREIPn7IuBvxVNC3DN1wCoPDTqhvsCUARfq6LV3kShGhI7WJLcCN1mTxvQr/M1M41NR7DTU8/M0KAGgFurBDnKCDIKIJ2CvHwZS5aRPqwApCma9K1dH4klqYHEDomF665piIEGalMDk/OHEkr9WlVFWC5xD7LRN9GVdW222kZ2/KjD3DsXiqO3UtVi9jT1C8KxdJQnZI1A9x5LhBC7L+ZrFRICDxD570Rjdl7FJVfLgsk324Z6aNm4TkT90qtnaVgsDm1Tp06KCwUTpxJKT+MnwtPGOIIOKrOzPo69utCmxCnmv5AG/oIjEJ9LqsSTJMvG5PclRCiVwqRikr6asnkFElxgVOXNymL0eQafqzqfsJJwIO9lEER95JzWH9PR3sxfN3VTZ/MfNh1JREAICtdOWty6meurd8uKdq76hAouYJbRURoSuVyhPoonauZ6GoHLRUahNKjJGYUYOHfmvM9ZuQbVoWEXv8US0T1/v7P9PaC2n8uJ+4r88WxPnGc28XQFkq3i/EdqqvtP7y1QgDk3qssGYMlgKVLl+Kjjz7CyZMnkZ6ejpycHN4fpXw0C1f4FRZJ5YL+AI1V0jQYUi6na7TmtAw/nX6i8TNDNHHFqjnuBMagqnUjhGgNE9fGKw0pRpjv7kWWfgsOak4C4l7yyzFFl2oyVX3i9GH5kEa89yuHKt4fnPoau00qkyPYywXxS3vj4aKegsfpUV9hCmSuGV2aOG6+uOwCqU6HZG5wQEVcA4kZhcjkCFhZhcX4+cwTrT6pqr6EADDm18solnF9fxT555gAEHstQqGQdpIGNlg2qr7QqhSXyHE3KafSBWhxp+TRD9tDLBZh+ZCG6NMgGN+/3Vhwn/ucsnJs1RmVe9iDRT2wf0o7DG4WClUYIdGSTahcDDan9uihyG7cuXNn3nZCCEQiEWQy6xi4pcL4n73IKoRL6Wvu5Regkt5BNZO0tjnsJ6Dp0Adtvjeq3E/JZa8FQDg61VclQujm82yBVuWDKb6uSwvn7CBGkVReqTQRcjnBtcRM1Az00JhDUAhV4fqShjqf+tAo1FvQjFeF89vro/hqHlG6qCkN+9flEze4WSjmlJpKcoqkgsJJkVTGOjDLKlgTByj8Txn2XH2BPdAekZ0ksBBhkhgzyAkwdvMV/Hs/FWc+6cQzz35/PE6r+RvQromnmB9uugshxm25gpMP0vBl//oY1lLYt9mamPn7DfwR+xxzSzX13q4OqFkqXNUL8cKP7zTR6ALAXbgLaeIARRBh/apegvsz9wZuaUZLxmAh7sSJE7obUcoMcwGdfpiGLqUJTbkXoLMjfxVtiKlLKCeXPnCdrrkCmiZuPs8WzHHHoJpGJKOAr03jPmT1IbtACi9XYWGl58ozWvd1d7JHkbS4XEJcep4E47bEonPdAEzqGFXm4xiLw3dSMGmbIh/azfnd9BbkNEUImyp2SZ/FAXMdMH5zzO+kaV87sQg+bo7IyC9GoVQmmIYjt6gEzg52IITwVvoVpY3NLDDM1JldKMXLnCIEemo270hlcvxbakb6PfY5GlT1YlMCLT/6kCfECfmpTt95Hf0aq/vdUSyD01oqjsjkhC3YvvFcfKUQ4hg/5i/2K1wGhO5NmioVcc2gZVFManqWWGrwm8FCXIcOHUzRD0opr9X0Y18zSQy5qmBV4YYQhQaBEbS0RfOMahuBS08z2Ju9LiQlMuRLZLwHplRG1EqXAIp0B4zWiyuU6ZpEQsKToUJc9xWn8d//OvNMCUKLNDdHO+QX8zXFzM2hPCkWhv18EfdTchH7LBMDm1ZT05ZWNL+X+o8BQIP5R/B0SS+9Akc0+ZkZW4hbPrghfjr9BF+91UBnW67WtkQm16mJA5T55ZKyCgUFM4XJ0Unt2qsoIa4sQQTPMwt4QpxIxJ9b3L6LoPlBBGiek5n5xTxNKcVy0LbgOXgrmX1d2cypDIbcggI480SpidP/CJoWvRn5xWqBdJaAwT5xqubSixcv4vTp05BKrUP1aOlw82QxmaN5mjiBnG0FHNv9U4EyRhGlGjgnezusHtaE3c74PGmi36rzaPHlMV6Zr8JiYXM5V+hS84vTglQmR0I6Pz+XIX5+gDJhKq+ckZzwfI8A4L02EWr7MjfH8mjiuD4YeRVch5UQgoT0At7N+3I8P9HtTE4yTG04CAjnQNl84rQxoEk1HJ7eXqugwRDEWVUXlegnxDFO/ScfpAkmC91w9ikA9bJVFWVOvf3CcN/h4hLt5YNinyl/86ISmZqgpk/SVOoXarlou965VXsSMypn0KEmIezx4l5q26QlfF9RwDAhEADealJNbZsuv0RzobcQl5ycjHbt2sHJyQkdOnRAZmYm+vTpg9atW6Njx46oX78+kpOTdR+IohWRSMT6xTF+KtwLUChRboFEKfQM+ek/rcfnJmoN8OSvKlRXcfeSc1AiJ4jlCAVnHwmHWnN35fqhaco2zyCVydWiVZmoSG37qJ+fYA+npuCRuylovPAor42QSt7OjtHEGecBli8xjU+oTE5w4XG6WuLidaefoP03JxA55yAepeahoLhErY2+tRYrShNnCNzr9VWuRGdgAwBEhygWJ7lFJYK/a07pvFJVvj54qT11QXlgkngD4AUk6IuqpthBZfxLOD5T6XnFakLc2tPKChCarnRLNBVRFGgLVOEuQF2d9LdgWBOa7kFC1yyjGOjzwxnsv5msdX9NCFmbDFUuVBR6C3GzZs0CIQR//vkngoOD0adPH+Tk5CAxMRHPnj1DYGAgvvzyS1P21WZgJqxQcXChFcmVZ9od0Ln7cF+rBhRwSyBpUstnFWry51G2n7vvNvt64tZYrX2TyoiaUKaq7XuRVciLFBKaTOcfp/O+r2P3+CbjHvWCUDtIPTSd+a6NFdjw3sZLZd53+8UEvPfLJaw//UTtO9h47ineXv8f6s/7BzsuJbCfcxMZd1l+CvP23Snz+S1RD8O9Xvdce6GXJu6NBoq8iUUlMkENFJPHTpdgZEmozhHV5MDc/JFCvkI/n3mq8xwJGaavWEFRIpcTnHv0Co/TdC8etAkh3GoCWQVSvbSu1oYhUzMjvxjrzzzhabwNndlCc+h+imVm39BbiDt27BiWLVuGN954A6tXr8aFCxcwb948VK1aFaGhoViwYAEOHTpkyr7aLLpWEbrSc2jaXTXfzqbz8exrrlDDNTtpMqdyZT5uOSBuhM+xGe3V9muz9DjmqggeXCHtTFwa2i79F/1Xn2e3FQn04XlmgdbSKMsGNxSsi8cIA8ZKscDNz0WI4iY9fMNFrOPUwhTiakIm/vfnLZx6mIYvD95D3bmH0fiLI/hw53XkFkl5EY2z99xC3bmHETH7gNqYf9cjsbEmjFmL05hEBSjShkhKZPr5xJUKNA9TclEs8LteLS3JpaqJUzWvWhKq5lTVaFJuourbL7LVBHLuokzTAq3fqnPl6yQF5x+9wi09o+2P3XuJYT9fROdlpwSTznLRVpZQ1Q86VkPdYGvGEJeO9PxiHL6Twt/fQFWckBB3+HaKQEvzo7cQl5mZiapVFdFLPj4+cHV1RXi4MgqmRo0a1JxqJJjcWAyqF/DvE1pjfPvq6NsoBEDZ89nY24kwqWMN9j13tS/jBQkoXx+5oyzezYV7G5Fo8ImLCvBQSy+iWvsVUDyEGBjz7b1kZQ6kLIF0COvPPNXqMO7mZC/44C+LT5wu36lVJx4hYvYBRM45iGE/X8SZuFdYcui+WjoYLneT1Fd5mQVS/HntBWLmH2GjM8sKk/xYG5piO8pa19ZYdK+niNI+dvclG02qyfQLKBcU2YVSnr8QQ9VSgUdVaLPkrOzPddR15foB3XqRrSaocS9vbVe6qh8pRX9Ssovwzs8X0XfVWb3ab7+UwL7mLqBVOXQrGXe1lDNUlcmn77iu1/krmmfp+Zi9+yYu6VHvWhVDleRpKvWGDb2FCQlxmp5r5kZvIS4gIIAnpH3wwQe8gveZmZlwcxPO0k4xDBfVyEyVC7B5hA/m9KqLKqUF4uJfafc707SICfNxxSc96uDj7rUBANsuJrDCDPeBbqjPhURAMGPQJ+p0zzWlNokrwDK+REKBE49S8wQf2NroFh2o1MTpGZ2alitB8y+PYfZuzcEC3/wjXKtVW+JhU+epW3dKczJnBk2RzeYuj8zkRHvMMZdrK+1VtYpCSNMU0ZcnUSwCVL/zjPxik+VLK++vu+iA9jxhunw6GUFYF6r+lBT9eZGlELTlRHeU6J/XnrNpQQBh6wLDxNJ0QZpQvY71TXBe0Sw9dB87Lidi7OYrBu9b3hrthu7tKOB7rskKZW70FuIaNWqECxcusO+XLl3KE+LOnj2LBg10pwyg6EY1e7umC5BZGei68aru/+voFujXKASfdK8DQJkYFwCuJypU8Zo0cZq0ftybVqFUxtauVMXDWXdWG67mJyOfG3mlOKam7PaHNKi7J3dSahu5CY8HNq1msCau+ZfHkFkgxY7LiVh5LA4Rsw/otR+gEACXH3mAt3/6DyceKHz2CopLEDH7AOb9VXZfNn3QR8DV9Nwxt6fYAIH8ZdoWAx6lKQJyOfPi1Mcd0afUV+6/JxnIk5QI/uaWqonyc9ee+kM1bZDqyCI4Ue/M7zy0uTJbfZXSSGFryVJf0ZTI5Jix8zq2/veMt73vqnOImH0AEbMPYOJWpbCllr5GJsej1Dz2Pvnhzhu8zzVpebRV8mCwVDcIVRhtoq6FkpAArK8M16VugPAHBgqBEgFry2MBH3VLQG8hbt++fZg2bZrGz1u0aIGVK1capVO2Tqfa/AtR0yqEqWlqSE41AOhQyx8rhjZmUzwM5IRT55RGlnJvQtybRE6hsMCoOu1OxylWmU4qkadCZU5U4a54dnF8wZgkwRn5hmlLuOHi35eWfwIUwuK9ZIWZUh9T2loVv7bvjj00qB8zf7+B7/99hAtP0jFq42UQQvTO2VdeuL/PtYRMpOaoRwBrehiY25zarV6Q7kYcmGuOGyUd5uOKGpySXM/S89mFiqOdmBXuTRaBVsbn7N8ftAOg/hsIpUDQdr6X2eq/t5ijqWTuIbdeZOPm8yxM3nYVs/64ie+OPsSOSwk4+SAVD1JykV0grbS5yLRx+E4K9lx7gc/2KoO2ikvkuJGYxb5P5Zjw1p1+wvueFh24hy7LT+HgrRS2NBqXs49eIV1loSUpkeELLfVxGawljkGo3JsQQuMp0FML1jLSV3C7oXewVgLHMdTSU1EYp3o6gObNm6N+/frGOpxGVq9ejcjISDg7O6Np06Y4c0Z7Rv5Tp06hadOmcHZ2RvXq1bF27VqT97G86JM/C1Bmpta1WtOliq7JCXA4elfh88aNcOIKdHeTcwSjn5j7FZOTjvF1UxUwR7eL1NoXQHEzFFJdR/q5l36ueCCp+tdpgltPs0WkUnssFitNtNp8UgAgt0iKpYfu63U+fSGkbIlfy8LpB2kghODvG0nov/o83vn5olobTQ8Dc5tTDU19oZqyBlDMgQkdlBrZvKIS9joWi5U36F91XAcVSctIH7iUVmhRvU4YzVmz8Cp6HevXC8/UtENcc3NyqZBXJJVj47l4HLiVjJ1XErHyeBxm77mFkRsvo/uK02j4xRFEz/0Hr397Eu+s/w8zdl3H14fvY8uFeBy9+xK3X2TjVZ6k0kVICi1etblgfPPPA0TOOYh91xULT+b+suzIA0TP/Udwn5m/87Vzk7ddxbaLCYJtuViLUK2pwoIqQhryHD3dHN5pGcZLmM/A9bPWB1XlA6CfFckclKlXL168wLlz55Camgq5yoU8depUo3RMiJ07d2L69OlYvXo12rZti3Xr1qFnz564e/cuwsLC1No/ffoUvXr1wtixY7F161acO3cOkyZNgr+/P9566y2T9bO8qPrEaXqEMRearuS6hqwgtl9MwOL+MTytjOqkKpTK1HySmBuJv4cT4tMLkFVataGsN5hhP/+HPZPa8rZ9dfg+3mpSlX34BHg6I12H+WuIiuaPKxAY4meRYQIz2/Sd1/HXjSSjH1eIXEkJJm69ykZtPUrNw/lHr9AmSnnDY35zR3sxOtbyx5G7TBCLuQ2qfFR/U1U0lafipuFYdvQhvi6tGGEvFkMkkoMQvmuBufF1d2QffJoir5uGV8EVTqJfQHGNM/6NgZ5ObPR6nqQEnpwH0ZDmodh84Rl61g+CnViE/TeTsfL4Q7Z9h1r+CPF2QUp2IVJyJEjJLkRmgRSFUhmevMrnJQFXxcFOhEBPZwR7OSPIywXBXs6c94r//u5OaqlSrIG31pyHnUikl//ktB3XMY0TaKDqcM/lxAN+aS3VNEmaeJGp7gOXmlPEq1xgCQgtroQQsgho8/lkrnEnezHcnOwxul2kmmUlSUATrQ0nB2VfP+paC8uOPrTYwAaDhbiNGzdiwoQJcHR0hK+vr1oOMlMKccuXL8eYMWPw/vvvAwBWrFiBf/75B2vWrMGSJUvU2q9duxZhYWFYsWIFAKBu3bq4cuUKvv32W41CnEQigUSinGg5ORWfG4YxkzJokjWYSXHsXqrWmqb6qKIbVPPCzefZqFmazoG7j6qgJCjElf5n8q7Fl+ac8/NwYk20hsCkgeCSlivBs3RlKpEwHxf8r1cdDN8gnJtt3fCmaBvFX5VxvyNtEY6qPBe4UZaXihLgGFTD7mfvuYXTn3Ri3zMCd4iXM880ZG5NHACM71CdDc54t1W41rbuTvboXi8Q/2iIpAaAm8+z2AeDWAR88WY9fL7vjskCTPQ96ndDGuLOixzcS8nBhA41WCFHY4JgEVDd342XI1EmJ6xmvH6IF9JyUyEngEQqA+HM22AvF9xf2ANO9mK2qgc3XVH3ekF4pyV/cVwklSEluwjJ2UVIySlESrZCuFO8V2x/lSeBVEbwPLOwdN4Ip7wQi4AAD2cEejkj2FMp3AV5OSPI0xnBXi4I9HKCk0CVGlNyLzkHcal5WLT/Lm8eMMQ+Ex6PPuSaIHBkj0otagCYu+8O1g5vavRzlQd3LcFIXAz18Vs5tDF+PvMUk0p9n90c1c9jaDUWbunEl6WWn+ISuV61wysag4W4uXPnYu7cuZgzZw7EWrJIG5vi4mLExsZi9uzZvO3dunXD+fPnBfe5cOECunXrxtvWvXt3bNiwAVKpFA4O6mbLJUuWYMGCBcbreBnwU6nPpumaYfJnAYp0FD4azIv6lMF6t1U4PvnjJqupWn9GczRjUlYhL0AAAPuUYsw+zITV1w8CAN5sGMITbISczLkO6fZiMV6r6Q8vFwfeqjjEyxmTX49Cdw2+VK2r++LJqzw01dMUBQBLDmmPDrRGVJO7MvKLWCRC93pBuF7q72MJOXBbV/dlhThN5cG4VKviKri9V0wQDt5KQZFUjtxSrZudWASnUu23tlyDxqB1dV9ceJKu8XNvV0d81ieafZ/KeYAUFstYbSK3l0kC0YiM8CESKTT7+cUywcUc4+7A1TwwCE1dZwc7RPi58QIlVJHK5EjN5Qh3rNCneJ2SXYSXOUUokRPFtpwi3NB4NIXbBCPYKQU9vnZPW7SyIRQWy9Dr+zNlKpxuSaTmFmHzhXjU8HdXW8iaC261oeIS9Uo9DIYupFpV90Wr6kofNm8BdyRD84CG+ijvH9x9EzIKEO5rWVk4DL7yCwoKMHTo0AoV4ADg1atXkMlkCAzkh8oHBgYiJUU4KjElJUWwfUlJCV69eoXgYPXkr3PmzMGMGTPY9zk5OQgN1e2Mb0w8nOzhaCdmV9+aEh02qObNvi6vQzYjbDF1QDVFegIKUxz33IDyoVIvxBNXnmWyD0h9/QiquDpgYd/6PCFu0rarCPZyZv11ACapr+J7sS+9KWx4rxkGrlVGTp+f01nruba93xIlcgJHezGrsdHld1WWepdcPJztUVAsw/TONbHsqGEBEWXB0U6MM7M6oeXi43rvw/gxiUR804exa6eWBVfO6lpopa2KWpqeUmoHeuLgLcW1PaG0moidWMRe/6byUSQcU7U2QlWET65P6W+XEtR8SkUQoVPtALX5qnSSFyG/VHjLKCjmPZy4v6qQlqSsAS0OdmJU9XZh8/EJIZMTpOdJeMKdQuBTCH4vS7V6khI50vOLkZ5fjDsCuRQZPJztWeEuyNOJFfIYoS/Y0wWeLvY6tSi5RVKzCnBbLsRjeOuIch/nakIWa82IX9q73MczBtxrf8flBIzQMM7yKsOFfr8RbbRr77VhzxE+Tz9Mw/DWVi7EjRkzBr///ruaRqyiUJ2EutSbQu2FtjM4OTnByclJ8LOKQiQSwcfNkS3sru2+w2ihypvDhtHqMTdzTTm2AGHNHvO9MukdGE1A8wgftYLsQgxtEaYW0HHhSToCVLSSs3bfwsTSBMWM6bZZhA+GtwrHFpXwf02IxSI4lo6vbZQf/rnzEi0ifDS216csji5uze/Ovq4IIW7jqOYI9HTm+UQJUSKTsyY7riaOJ8SZX4ZDg2pe6FjbH1VcHXmCiCZ8OSk5vnorhn1dI0B5A2a+Fyd7O/ZGrY/Wujzocu7matcBwNNZOSfSOL6t3AeV0IKLEdxEIsDNUaGJK5ERpOcLm8kHNq2Gn07zte+mrKVqJxYhwNMZAZ7OaKihDSEEWQVSrabblOwi5ElKkFtUgtyiPDanoBDODmIEe7mUmmpLzbYc022Ql7PZq3Z8vu8OBjcPrXATckXArf/6UiA6nqG8QTHuHMVBrUB3rH23KS+4zVC8XBzQoZY/Tj1ME0xOb24MFuKWLFmCPn364PDhw4iJiVEzSS5fvtxonePi5+cHOzs7Na1bamqqmraNISgoSLC9vb09fH2FQ5EtBX2TjjI32ld5ErUHgKujnd6h2f6lwlKepARpuRKtq3ChvHTMtPN0UVxSl55msMENADCyTQT7elaPOvjqMD/Sk8k4/8/09ui+4jS7XdWhtUgq45hTlX384PUoJGUVYkhzw7SmjFlYW6TZZ3/e1vgZw0/Dm2LcFu11YhnaRvni3CPNJjVNuDvZ652Mlfn5dD2Toj49hI2jmuOPK89x4JYimXdydhEv2aW5U4wACo3UplEt9G7PTTXA7X/vmGB8gGu8thn5xZwAAhNp4kr/lyWnF+MPKBUQMDX9NHuuKlPzhPu64W5yDgqlMszZfUuwvZDm0pRCnD6IRCJUcXNEFTdHRId4amyXWyRlNXeMYKeq3csskKJIKsfTV/l4qis5upmp/dlhfDPQeDlXlx95gEmdogxORWVMCotl7P1F8V7zPCtv3juuBlgsEpVZgPusd138Efsc77WJwHdH4wBYVuATg8FC3OLFi/HPP/+gdm1Fln9NxdWNjaOjI5o2bYqjR4+if//+7PajR4+ib9++gvu0bt0af//9N2/bkSNH0KxZM0F/OEuid4Ng/KFHHUwmalLownd20F+Ic+as/LZciEf7Wn747VKiYNtFB+5BLBLxTDvM6ZtzNFrcmyX30pjQoTreaBgMRzsxWpSa+xhzsKsj/0bDjG9ypxpYdeIxJCVy5JcKMlw1d6CnMzaMbK7XWLkwgqA2X6gCPUzV2vxOOtfh5/1bPrgRtv33DDUC3FG/qhdcHe1AiMLc4GQvxvnH6RgvIBDun9IOXi4OeJZRoLXOpZ+7I5qEKfz99Lkdjtp4mfc+T1LCX3GaX4YzmEiOzxZ3bgjdo/7Xuy5rTr2akAW5nPByqBkTbUJ4kIZoQl2m3nohnmqmRmbtI4IyKvfE/VS2jJ3iM+UYhcyp5hbi9MXD2QEezg6ICvDQ2IYJyNBmuk3Lk1iML9zHf2iuCGMo3//7CD5ujhjZVnd6J1OhquXdevEZ5r4RLdhWRpRuHd4uDrz624ZSnt/z/deq4/3XqgNQ+vMZmqqkIjBYiFu+fDl++eUXjBw50gTd0c6MGTMwfPhwNGvWDK1bt8ZPP/2EhIQETJgwAYDCn+3FixfYvHkzAGDChAn48ccfMWPGDIwdOxYXLlzAhg0b8Ntvv1V43w2FuzLWJhwzN3Ch8Odm4VVw5O5LNp+cNrgq6O//faSz/Rf77/KEOOZB6efuhKgAdzxKzUORVM4KEdwHhkgkUnM8Z5xHueYjLlyfqNUnFUl39c07pA3mGNqcaXP1WH3ZC5RpAYBfRjZD6+p8AS/Q0xkzutXWeKzu9YIQv7Q3iqQyzP/rDnZcVgjTwd7OcLK3QxUd+fEu/q8L+wAuaxkpbpCMdTzK+XDzPOkykTrZi1EzULlav5aYZVDQiyFwTUUezvZsQmJvVwes4CSi5sJqCTn7ckukCTlyc2EEtMdpeegaHYj9N0s1IpwfVuiaMiR629LRNyBj7OYrvHJYlYX5f981qxCnGihXXLoYFwpIIRy3jj4NQvR2kxGiUah3mfflwrgr+XuY19VKCIOfgk5OTmjbtq3uhiZgyJAhWLFiBb744gs0atQIp0+fxsGDBxEernBaTE5ORkKCMjliZGQkDh48iJMnT6JRo0ZYuHAhvv/+e4vOEcfA1Uhpu5VqyxU37816+KBTFNbpEWpuJxZp9QvTRE6RFEVSmVJYEyn7JCmR6b0SYoRAL1cHDGqqno1eKL+SMTQFdqwmTvODPj1Pd444J3s7zHsjGpF+bjgwtR3uL+yB+KW98XqdQF5+MkNwdrDD0rca4PHiXoj7sqegn8zEjjVw/KMOvG3c70VVs6kv+i4iLBWuJk01Pcyqd5rw3jvaiXmLClOYTJh5wE0Vwv09V7/ThBdhx4UV4oTMqQDmvVFP43lFIuCt0vlUIiO860F14dI8gi+4mkobaak42In1ToNhLKr7V5yT/PnHuqvSmAqhRbImrTTT1k4kQpNw7zKd79iM9pjWuSY+7VO3TPurUjtIoeU1tc9sWTBYiJs2bRp++OEHU/RFLyZNmoT4+HhIJBLExsaiffv27GebNm3CyZMnee07dOiAq1evQiKR4OnTp6zWztLhap60PUMZB/RDHH8DhqreLpjZvbZaJKkmDM1InVVQjFaLj6PO54fZPDwikYgV4g7cVPZJ0xh6xyjMqv2bKOtjfjOoIc+/z04sYitBcNGk/TIE5hjaNHG6tFkb3msGABjVNhInZnZEvRAvo/qf2IlFGrWONfzdUcPfHdc+74qm4VWwfDDfTVw1MERf7Hk+cWU6hMXwuoo520HlumG+W6b6gcSE9UO5vqGTOiorSGgTlJn+/vdU6UfJXRzVCvRA12hhv2ARRKx/o1Qm512Xqn6g0cF8v7OKqiZiSegTNGNMvh/auMLO9eHO6xV2LlWEsidoWoTLOebUfo2q4qu3YnDkw/aCbTURFeCBD7vW0mjZMRSlYsLy5oTBy45Lly7h33//xf79+1GvXj0137I9e/YYrXO2jDMnb5O2FA+Fpb5LiUZIRvvB61E4rlLL08FOxPqLOdmLeRfxd0cfsj53XB8cZoWVnl8Mn9IoQU0j+PGdxiiSytW0VU3DquBRqiLSbG6faN73wVBTiw+MvjAPcE0ZwZ9oiUztFROEH99uYhaNxc5xrXA5PoMtDl/FzRG7J7ZRa1eeNBEM1irDXfu8KxIyCtBQxaTSTqUsDyOwMtfgg5Q89DByBUHG/Bnh54Yzn3RCTpEU1bxd8cV+RW1MbZVNmAVdYkYh0nIlPJMO8/OqCqZcuEEbXE2Tv0quR2fOHPT3cELbGpaRX6wiGd02EmtOPhb8zMXBTmMqp/sLeyApqxDnHqfj8726A6EYHO3FiKnqhVsV4Gv1MkdiUn9Pbfi6OaoljS/S8F0yaws7sQgikQhDmqtXY6poGK25ttRb5sJgTZy3tzcGDBiADh06wM/PD15eXrw/inHQ1wQ3vr3C8dLQjNRCNA6rgmMz+KY5bq3RgSpmzl8vqPsqiEUiTOoYBaC0pisj3Gm4b4hEIsGxzn+zHjaPboHbC7rjvTYRvAzamvpTFpjVoKboVE2RbNM618TqYU3NZnJqWd0XH7xeU+f5vx0knMDh5vxugtsZylqezJKo4uaoJsABCqHo/OzX2fdM4ABjdmXqlZqKUB9X1Avx4iXYlWiZvz3qK5NWayrbpElYF4mUQlwxJ3hndNtItd81kDPHjn3YQafvZWVEm8+TtqhJZwc7VPd3x/BW4TrnFpcQbxf8OUl98WUqRvwiXN3G1NSvqpANuFrx+FcFgm2Z79kSouIZqlZRRrzqqlVe0ZSp7BbF9ERwskLna4lo83bh52UrL25OfIGKm1hVH/OKSKQUQJ+lF5QplQagOEb7Wv7s+461/TH19Si9gi4MwaE0d5GmjN6qCrpfR7dAmxq+RgmqqAjqV/XCZ73r4q8bSXizYQgWHVBUnnDVYe7lanYs6F5qNLjmc8aE1qq6D56+yjdpLijuV6lv8IWPmyPCfV3xLL1AoyZIHyFOKpNrjVYe3joc8en5CPZyUcvZSAHebhHGFrLnIpTbL35pbxy6lYwlh+6rVUYBFOXV2tbwYzWjpz/uhPbfnDBJv7mcffQKL3OKNNYXNhWMYNanQTCuJ2YhI78YxTLha1nGCnEV1j2dcAOdpDKCCnad1IoFdYXC5TWOyUdbNQbGJy7ZwAK/mlB9GIRwcu7oI5CJoBT8uA7l5dXmiEQizOhWG34eTpi77065jsWFMaUlZxehSCrj+QwVSWW8fFsArEqAY+CGyrs62qOGv5vOoBBuYk5LqNhgbKpVccHbLcIgKZGxD2Hmt7+XbPx6yYkZ6u4OIpEI9UI88SQtn00JowkmBdCLrELeA4X5bbT9no6lZcoepeahTQ3f0nOrt3OwE+OLvka2I1shY9pFYsPZpwAUqY261wvC0bsvMblTFD7sWostf5YvKUFxiVyjH13PmGD0jAnGzedZIAToy0kLVMPfnVegPkzA59dU9Fp5BrGfd62w8wH8iNMof3dcys/QuFhiXAssKbCGe8+XlsgBCwpSpUKchcIVejT5DgCAZ6kmztCgBE1IVCYWN3FinqQEPesHafcLEPGz4nM2G4V+jati1YlHaBflr7uxHnBXpJfjM/BaTeVxt1x4xhtrVIC71QlwqqgWNNcE35xqqt6YD5FIhCUDYnjbmMXHodspagJ9eeAWTFddzPw5qS2KSmQ6HbCfZSjM+ozbhKoP3ZDmofjz2gt0rhOAp+n5eJKmaC+CCNX9lJqiM3Hmi1C0Fj7vE40JHWrA08We9YVigsOcHezgVXrPZf7rgtn38qdd8PXh+6haxQUxVc3nepSeX4wXWYVay6IZG0bzLRaLWDeCLw/cQ68Y9dKXjMHHklLc2IlFEIsUlhlLC/ix7ieSjaBNcGBuJPmSEq3O0foS4s1Xs3Mf5pISGVYPa4K/P2incX8RRAj2cuHlGSv9wCh4Ojvg/OzOWDZY2NfLUHzcHNm+qiZGXnea7+DM9Q+s7HCvuZwy5pqzNgZwIqQzC3SnldGXawmay8452ov1iqBj0o9oio5rVd0XZz7phDXvNuWnlREpfAOZXJHlLc9nK/h7OBm99JW/hxO+GdQQ07vUMruf6Yc7rlfo+eQcEynjd6jpvqKMTrUcIQ4QztdoCVAhzoIZ3746vF0deI7NqjDJEuXEOOHP9nZivNEwhH3P1WgXSeUQiUSIqeaFf1VykzEw865jLb6m7L/HZfONE8LYmeRrByqiXJnvL7tAim0Xn+GVSn64p2mWXa7HmHBTjGirRlGZqO7vDo/S+WRMvzhtdYj1hZnnquluuM+5UB9XONqLeaZwhjqlea4YLYJlPR4pADSmieFiLO3ZpfgMpJS64Dx9lW+UwDhtcM2pH5UmOtcUzMPmibMw6cRRS75Gc2JhXxOFy+yedXDt8648vzRVuA7q8enGETJ8OE7NmQVSNp/Vd0OU2q/q/u74/u3GannIGJ86J5WUIHdN4GdkLJi+MvnBGn5xBJ8K1Et1sK9c00UocpPBg+O5u/9mUgX0xjJg0mxwy1OVF3sjpGthgiCWHlLUHNamC3AUOJ8yQtWyHkAUJWuGNRHczr3HquaBLA+tlhzH4dsp6PTtSSw6cA+EEIzbfAURsw8gYvYBRP3vIC/XZ3ngauKcOQnqhYrdcwU+S4K5/6flGSeI0Fjo7UjFlLLSxYgRI8rcGQoffdTJXOfPW8+zUa2KC55nFmLtu7qrNGhibPvqbPoQqUyOT3rUwdTONdV8hN5sGII3G4Zg9u6bbGkopjeq2gBNedgsAcZpvEjHCsutjNUPLJVWkT64kZgl+Bm3HI61+wEaAmPiOfUgFcNbhRvlmFzNsayMLg8xVb2w5+oL+Lnz3RSE7hAO9upbmQcQq4mzrOcjBXxhn0t1fzc2+0CQHiUUDYFJALzpfDzebhGGI3dfsp+VyAkmb7+KpuGdNZ6XEKLXc0rGSRvCfY4sPngPn/WJ1tjWkmDyn+65+pxXI9zc6C3EjRw5Eu7u7rC3t9foeyUSiagQZwZaRPrg0tMMXtqE8kx2bgki5qavzcmbm+eNmXfd6gXyat5ZSmFpIbiauCvxGRrbjeHUiq0MqP4kw1uFI8TbBT3rB/HSX0zkVBao7Ex5PQrfHnlo1OuVK8RpSmWji9frBGDB33f18mnjJWounZCMds4SywZRtDO8VQSiAtzxep0Aows23MwHmhz2k7ML1Z4nf157jg933gAAbBzVHJ1qBwjtyiLnaNe4Jf1+PvsUtYI8MLhZKKctEwSh/zgqgvohnriakKUW/Gdu9P6a6tatC0dHR4wYMQKnTp1CZmam2l9GhuYHIMV0MKtzqUxutIePb6mzfw1/dx0tFRo5Biblgb6lviwBRhMnKZFj4NoLGts1s6DVlzFQXYyNahuBiR1rIMLPjacV4OYsrOwwrgvGNDtyc+6V1feIefDlF8tQUFyidZ5zAyWU5lTFKwtWiFMgrO13dbLDon4xeL1OoEnN4ZrcceLT8/HZ3lt4XFq9ZsWxh6wABwCjNl7WeWzCEczEYhGmvh7FfvbJHzdxJ0lZsYIxsVqaJq5/aXUcbSm/zIHeQtydO3dw4MABFBYWon379mjWrBnWrFmDnBzL9XWyFRjTZXEZV/lCfP92Y0zvUlOvqgi8MkClV5SLgOnVUpGWVmvQVG4H4Ne8rCx4qERFqppF1o9ohm8GNqjwepLmhJsY11hwH0Y5RZoTd2vDk5POYuO5eKTmluaFFHjQjSut4sJFdTiWFvlHUbBrQmve+061/fEaJ7DIlL/aB9uvCW7/cOcNbP0vARO2xOJMXBpWHIsz+NgyFcFM1XT8zT8P2NfMQsOSUowASmtUXKrmUozmwCCFZcuWLbFu3TokJydj6tSp2LVrF4KDgzFs2DBIJJbl7GdLmOLB0zbKD9O71NLop8FFKFpUtZajak4uS+LoHYUfSJ6GyhgHp76GeW9EC35mzYxWMQ+raua6RgdiEMfMYQso55JpVFZlLdnj7GDHasfTciVsOpwCgWvWVcC9oVEoLYloDdQL8UL3eooo1WAvZ2wc1YJ3D+YuCJio+ooiLjUPwzeUrWyXXEewAldDzAh8FibDsdpvmYWps8tkdXZxccGIESOwYMECtGjRAjt27EBBgXAdNIrpYTKymyv0meuDQ0q7wF3ph3g58xzlLQ1dqVmiQzwrpebCXeU3saxbk3lg55KJzFbleQCMahsBQJH8m/nt3AWSfHN9VJnzqdZBrXxXc+Vh9bCm+G5IQ2we3ULtM+7vvXdyW5P14cniXjg5s6Pe7QetPa/1c6ISrNClLj+dSmpOERIzCpCaU8S2NXYqqfLiW+q2ZKwk4MbCYCHuxYsXWLx4MWrWrImhQ4eiefPmuHPnDqpU0V42hmI6GCHqwctcs5zfz90JfRuFoG+jEMGaiwUW5kOgiouWqNO3W1RuTdSBqcrEzZYcfFJRMK4JN59na62UYgjcwunlidJmHh47LieyfQvxUk8/5OumdG9gNHaq7g0Uy8VOLEL/xtVQU0DT5ufuhMX9Y7BsUEM4O5jO818sFiHcgFJgl+M1J7QG1IMVokM8MaNrLfbz+ym5eO3rE2ix+DgyCxQR4pbmE8e4NBjrvmAs9L4Kdu3ahZ49e6JmzZq4fPkyli1bhsTERHz99deoU6eOKftI0QETcbafk9Onoi//lUMbY+XQxoKfqVZCsDRqB2k2S/RtVFXjZ5WBukGe7Gtf1SobNgg3Cu/m82wtLfVHzlHqlSc1AbdE3IkHaQCUee24cF0ZGPNtddUAJct6PlIM4J2WYXiraTW9rQNNwrz5VTz0xJjWB1WfOACY2rmmYNtlRx+otbUEmIXQ01eWlfRdbxvX0KFDERYWhg8//BCBgYGIj4/HqlWr1NpNnTrVqB2k6KZ7/SA2T5sl1XWLDvbE3eQctLPwjP/LBzdEu69OqG2vX9WTV2y8MiIWi7B/SjsUSWVqJjdbpBZH+2GsdByM7i3c1xWRfmWP9O0dE4wpv/Gdz4VcKPh1lxWfl+e8FOvmf73qYqQeEaRC1PB3w2MjVKphjqEqmL3RMAR/3+AnE08vrZRjaSlGqlVRar1vPs+ymAwMegtxYWFhEIlE2L59u8Y2IpGICnFmoFWkL/vakpLqrhzaCLuvvsA7LfQrum4uuHnxuPz9QbtK6QunSn0zFuO2RBpU88LN59lGWxAxpiR90vVoQywWoXGYN64lZLHbAjydBNu6OtqhoFjGE95CvJyRVFpqiWI7NIvw0Rkk4OPmiAhfV1zlXFsAMKBJNV7kqDYO3UpGT4GC9o9SlW4+qnFyjgKBcwXFimAdS4tO9XZVLnJf5lhOIKfeQlx8fLwJu0EpD1zziSXVdasZ6IHZPa3D1P5x99q8m9WKIY1sQoCjqMMtUZWaW4QWXx6Ho70Y1z7vWqYAHcZR2xhXk2r1DE3VNHaNb40Lj9PRr7HSHcCRk8BZRO2pNoWuIIGhzUNx8al6nlduahsA+LRXXTg72mHlsThsfb8FHqfmY/L2q4rP9t4WFOKuJyrdElTvqX4e6tp/JjLcEu+/TGJ9S7J4WW7IIEVvuBOU1kYsG6qRmsYub0OxHphFkVQmx5cH7gFQmFZPPEhFnwaG5ztk4hqM8VC6pPKgtdfwcK5f1UtNw5qWaznaA0rFou3K69+4KqZ3qYW31/+n9tmgptWw7b9nyJOU4MCU19jAtXdbKixz3ETgjAZNG6rm1JFtIrDu1BPBtpYWnQooNYeWJMTpbXW+ePEiDh06xNu2efNmREZGIiAgAOPGjaO54syESCSyyIvLmuDeMD7oFIWWkZWrOgNFfxjt1pTfrmHfdaW/Tlk9FZQ5ssrbM3X0yePI0LqG0u3CApUcFBOiuoD4sn999vWgZtXgaC8WFPScHexweHp7nJ31Oi/zAHM8bnm+Ig3lqPKKpOxrVRNpsEB0NYMFynCwL13gWVL5Or3vAPPnz8fNmzfZ97du3cKYMWPQpUsXzJ49G3///TeWLFlikk5SdEPL6pSP9jX9AQCezvaY2b22RaryKRUDk0JANeXKH7HPy3Q8xifOFJeUJk2cEJ7O6ul/KLaB6lUytHkYvn+7MSZ3qoHW1RXCfVmuT9X7JFtJhANX+Av21t/CYWnRqYDpk4GXBb3NqdevX8fChQvZ9zt27EDLli2xfv16AEBoaCjmzZuH+fPnG72TFN042IsBC0/lYcmE+bri/OzX4eVCH3S2jqacV6cfppXpeMzt3hQPJXs7A4Q4em3bLKqXnlhUWgpRoO61oQxtHspmR3ialo8AD76gxtSmDvBwgp+7cCCOEJYoxDEWL31MxxWF3pq4zMxMBAYqsyyfOnUKPXr0YN83b94ciYmJxu0dRW+yCqS89xZ4/Vs8Id4uFl1ZglIxqPpHlhfVbPXlYXF/fvk6wzRxynHR24Ntsebdprz3gpaGMl4UjcO8tX7OlDOMCjAsOtsSfeLEpX06VcYFnSnQW4gLDAzE06dPAQDFxcW4evUqWrdWFuvNzc2Fg4PpVnqZmZkYPnw4vLy84OXlheHDhyMrK0vrPnv27EH37t3h5+cHkUiE69evm6x/5qZeiKfuRhQKRSeaauiWFTnj42CEZ1LfRvzACnsDkmm1rO6ruxGlUqJPkmljiEw3nmfx3sc+y8THfyjcsFI0pLdZObSRcH8sT4Zj3ZYsyWKj9x2gR48emD17Ns6cOYM5c+bA1dUVr732Gvv5zZs3UaNGDZN0EgDeeecdXL9+HYcPH8bhw4dx/fp1DB8+XOs++fn5aNu2LZYuXWqyflkKVVxpolYKxRh4GFsTV/rfGJo4riM5YJg5lVumyRIfkBTzUtZrgpsQffHB+7zPpnKSUz/RUOnAQ6D+L2CZmjhGGNYUxGEO9L5bLVq0CAMGDECHDh3g7u6OX3/9FY6OSsHhl19+Qbdu3UzSyXv37uHw4cP477//0LJlSwDA+vXr0bp1azx48AC1a9cW3I8R8gzJcSeRSHhRtjk5OWXveAViSSsDCsWa2TymBYb89J/RItCMqIiDvZ0Y9mIRm9TbEE0c9x5BA6AoqpTVJy4qQHPZQrkeBZk1FZS3RJ84pvTW88wCM/dEid53AH9/f5w5cwaZmZnIzMxE//79eZ///vvvmDdvntE7CAAXLlyAl5cXK8ABQKtWreDl5YXz588b9VxLlixhTbZeXl4IDbWOAui5RjYBUSi2SuOwKni4qKfadgcDtF5clD5x5eoWS8faikjqJmHe8PfQ31E8wteN1W7U0VIvmGI9qPpIlgdTyEz6CHGujsK6JEsU4piE2fdTcnW0rDgMrk7m5eUFOzt1ydnHx4enmTMmKSkpCAgIUNseEBCAlJQUo55rzpw5yM7OZv+sJVija13174dCoZSdLWNaoFGoN7aOUSwepTKCf++/NPg4hM0TZ5yH0s/vNceDRT2wZ1Jbg0xO9nZi3JrfDWdndULfRlV170CxeFR//u71AoUbAqgVqD2wwFgyU9xLpYDDLU/Vs36QYHs3R6U8wdUWW6A1FbVLFz+GBBSZGrOWmJ0/fz5EIpHWvytXrgAQjqYhhBg9n5eTkxM8PT15f9aAJpU0hUIpG6/V9MfeyW1RI0CZlf6Lv+8afBxWG2HEW5WTfdnmu6ujvcZawRTrQ/Xx5+2iWZEyuVOU9mOV4wId1TaCff3VYWX5wuhg5fPz8z7Rgvu6cnxQBzerxr62RE0c43teIifKgCUzY9Z8Ch988AGGDh2qtU1ERARu3ryJly/VV8BpaWm8tCe2jKO9WeVxCqXSwi3SXZYkn6bME0eh6Isuy2Z5Ls/hrcKx8Vw8ACA+XRnA4MB5LrlpMJtyNXFcP1RLrD7EDSSSyuVwEptfeWJWIc7Pzw9+fn4627Vu3RrZ2dm4dOkSWrRoAUBRBiw7Oxtt2rQxdTetAtVC2LTANYViHLjC14usQoP3lxvZJ45CYYgO5tfH1SaI6eOfVla4z59HqXkokcnx5FU+biRmKdvYC3eOm5fxGqd9TpFUoLV5UV3QWUJaUQvogm7q1q2LHj16YOzYsVi3bh0AYNy4cejTpw8vMrVOnTpYsmQJG3SRkZGBhIQEJCUp6h8+eKBQ8wYFBSEoSNg+b62oCnEUCsU42HFW32WJAldaU6kURzEuMdW0C3HdopWWKl0uN47leIa4OPKPPezni0jL49dS1xTAwK3/y+2+psop5oT7nC2xEE2h1Tz5t23bhpiYGHTr1g3dunVDgwYNsGXLFl6bBw8eIDs7m33/119/oXHjxujduzcAYOjQoWjcuDHWrl1boX2vCMoaOUehULTDrTnaKybY4P3Z6FSrudtSrIkPSn3dZnStBa4YdOWzLljLqdTQLToQXeoG4OPuwim5Pu8TjSBPZ42+a9rwc3fCoKZKf7aLTzPwJE1pVv1uSEO9jsN1M6tpYIWHisBOLGI16sUWIsRZhSYOUES/bt26VWsboqIuHjlyJEaOHGnCXlEoFFvg4+618c0/D/DbpQS80yJMTQOiDTZPHPWJo5iAj7rVwuBmoQj1ccGne2+z21XrlNrbifHze801HifCzw0X5rxe5uv0m0EN8Xvsc8HPXPQMvKviZvlJ6+3txCgukSMhvUCtTqw5oGvDSoIlXEwUSmWlsFjGvn7jx7MG7cv4IlERjmIKRCIRwnxdjbJIMNVCw1OHG8KG95qhdXVfLO5fH6E+LgCAltV1lwozB0zwRWaBZfjsUSGukqDqk0ChUIxHv8Zlz6tm7DxxFIomzH2FRfq5CW7XVbu1c91A/DauFapVccUPbzfB8FbhmNa5lim6WG5aVWdKb8l0tKwYqBBXSdBXXU2hUAwnSsU/5/vjcbz3mnJGvcqTYGVpWxqdSqns/PxeM8HthgTeNQr1xsJ+9Q2qRlKRMM/atFyJjpYVAxXiKgnerrR2KoVSUSw/+hBJpelGTjxIRfX/HcTr355UW51/dUhZEJz6xFFMTaCned1qgsx8/oqA0aifeJBq5p4ooEJcJcHZwY6nyqbPCwrFuCwZwK9TyQhsjKD25FU+6nx+mBXuAH7eKzonKaZm7GvV0a9RCNYMa2KW8wulMfFzt/xgBUNggi90+flVFFSIq0TUDqRFrSkUU/F2izDe+9eXnULE7ANqxbDbLP0XXx64i4z8YjxKzWO3U584iqlxcbTDiqGN0bMMqXCMgZ1YxAYmMByc9ppZ+mIqGod5A6B54igmwJ7miqNQLIL1Z57i68P3edvo7KTYAmc+eZ19XS/Es9JlTmD8+8pSgs8UUCGuEsHNuE0X/RSKedlxOZH3/tg99frPFEplxq4SRvM4skIc1cRRjAwtvUWhmJb7C3uUed/49AIj9oRCsVwGllZveP+16mbuifFhLF5MvjhzYzUVGyi64ZpTaZ1GCsX46Ko/SaFQgC/718eYdpGoE1T5/LQZZcnFpxlm7okCqrqpRDjaK39OObEMez2FUtmoF+IpuP3DLpaZnJRCqWic7O1QN9izUqbVsbScrFSIq0QM5BQgzpeUmLEnFErlZd3wplg/ohmeLumFFpHKTPSv1wlAiJdmJ+7+5aj6QKFQLIMWkT4Y2LQaetQLMndXAFBzaqUiOlipIZBYiL2eQqlsVKviimpVXAEAXpxcUTHVvGAv4Jf6doswhPm4YuxrkRXWRwqFYhqcHezw7aCG5u4GC9XEVSK4qmtLqetGoVRmVI1FH3SKUmsTU9ULEzvWEBTwKBQKpTzQu0olJczX1dxdoFAqPbN61oG7kz2mdq4JABjUrJqaIEf9UykUiqkQEULvMNrIycmBl5cXsrOz4ekp7NBsSdxJykZSVhG6RgeauysUik0gkxO1fFgFxSWInvsPAOD3Ca3RPMJHaFcKhUIRRF/ZgwpxOrA2IY5CoVgGWQXFeJUnQVRA5UuzQKFQTIu+sgcNbKBQKBQT4O3qCG/XylX8m0KhWBbUJ45CoVAoFArFCqFCHIVCoVAoFIoVQs2pOmBcBnNycszcEwqFQqFQKLYAI3PoClugQpwOcnNzAQChoaFm7gmFQqFQKBRbIjc3F15eXho/p9GpOpDL5UhKSoKHh4fJ6sDl5OQgNDQUiYmJlT4C1pbGCtjOeG1lnIBtjRWwrfHSsVZOrHGshBDk5uYiJCQEYrFmzzeqidOBWCxGtWrVdDc0Ap6enlZzgZUXWxorYDvjtZVxArY1VsC2xkvHWjmxtrFq08Ax0MAGCoVCoVAoFCuECnEUCoVCoVAoVggV4iwAJycnzJs3D05OTubuismxpbECtjNeWxknYFtjBWxrvHSslZPKPFYa2EChUCgUCoVihVBNHIVCoVAoFIoVQoU4CoVCoVAoFCuECnEUCoVCoVAoVggV4igUCoVCoVCsECrEUSgUCoVCoVghVIijUCgUCoVCsUKoEFcBXLlyBUVFRebuBoVSbmwhIxGdr5TKBJ2zlRsqxJmQJ0+eoG/fvmjRogV27dpl7u6YlMTERPzxxx+4evUqpFIpgMp988jIyMCrV68AAHK53My9MR3JyckYNGgQdu7cCaByj9WW5itA52xlhc5Z24IKcSaAEIJJkyahZs2aEIlE8PLygru7u7m7ZTLmzJmDWrVqYdmyZWjTpg0mTpyIJ0+eQCQSVcqHwqeffoo6dergp59+AgCIxZV3Gm3YsAG7d+/GihUrUFBQADs7u0r3ULC1+QrQOUvnrHVji3NWE5X3SjYTe/fuhZubG2JjY3H+/Hns3bsXdevWxaFDhwBUvpXuxYsXsW/fPvzxxx84ceIEfv75Z8TFxWH48OEAAJFIZOYeGo+srCyMGTMGx44dQ1hYGP777z9cvnwZQOX7XRnOnz+PIUOGwMnJCV9//bW5u2N0bG2+AnTO0jlr3djinNUGFeKMAPeiSUtLw9atW3Hx4kW0bNkShYWFqFGjBjIyMlBQUFCpbpCAYkLJZDL07t0bzs7OePfdd7F06VLcvHkT3333HQDrnlTcvru4uCA8PBxz5szBsmXL8OLFC/z555+QSqVWr8FQ7XtJSQkAIDg4GEOGDEGbNm2wa9cu3Lt3D2KxuNKM1dbmK0DnLJ2z1oetz1ltUCGunBQWFqK4uJh9P2bMGAwYMAAAIJPJ4OLiAj8/Pzx69Aiurq5WrdZmJhJ3DAEBAXBxcUFBQQG7rVWrVpg5cyYWLlwIiURitZNK9bd1dHTEtGnT0K9fP3To0AGdOnXC6dOncfToUTP2svyojpMQAnt7ewDA5cuXUatWLfTv3x9BQUFYu3YtiouLcffuXXN1t1zY0nwF6Jylc5bO2coOFeLKwZw5c9CuXTv06dMH33//PXJzcyEWi9mLiLkRdunSBfHx8UhISLBaX4zly5dj8eLFAPj+JJ6enrC3t8fx48fZbSKRCO+99x5cXV2tdmWv+tvm5ORAJBLB09OT/X2nTp0KQgj27t2LV69eWeXKXtM45XI5Xrx4ATc3N0RERKB58+Z44403sH37djg7O+Pff//l3VitAVuarwCds3TO0jlrC9jWaI1EcXExBg0ahL/++guffPIJQkJCsG7dOrz99tsAlDdM5r9MJoOvry8SExPN1ueycvnyZXTq1AkzZ87Enj17cOHCBQBgo9kGDRqE4uJiHD58GKmpqex+wcHB6Nq1Kx4+fAiZTGY1K3tNv+0777wDQHHTYG4iYWFhGDx4MK5evYr9+/ezn1vDQ0HXOMViMTw9PeHg4ACRSIQ///wTixYtglQqRUxMDKZMmQJHR0erHmtlnK8AnbN0ztI5a1MQisHcvXuX1KxZkxw5coTddvbsWeLi4kK+/vprIpfLCSGEyGQyQggh6enpxNHRkezfv5+33RpYuHAhGThwINm4cSPp1q0bef/999nPiouLCSGErFq1itSqVYv89NNPvH3btm1LxowZU6H9LS+G/rZFRUWkV69eZPDgweTmzZtk69atZNGiRWbpuyHoGichhBw/fpwEBweT+vXrE29vb/Ltt9+SdevWkUaNGpFVq1YRQqzjWral+UoInbOE0DlL56zlj9FYUCGuDMTGxhKRSETS09MJIYS9oJYsWUKqVKlCHj58yGuflZVF2rdvTz766KMK72tZYcb07Nkzcv78eUKIYnwtW7Yku3btIoQQIpVK2fbvvPMOadSoEVm3bh3JzMwksbGxpEmTJmTHjh0V3/lyYMhvy9wo9u7dS6pXr058fX2Jo6Mj+fbbbyu+4waibZze3t7kyZMnRCqVkujoaDJu3Djy9OlTQgghSUlJZPDgwaR9+/akqKjIXN03CFuYr4TQOUvnLJ2z1jZnjQEV4srAtWvXSL169cgPP/xACFFeYMXFxSQyMpK9kJgbZklJCalZsyaZMGECuxK2Rh4/fkz69etH+vXrRzIyMgghhEgkEvazuXPnEjs7O9K0aVPi4uJCxowZY3Xj1fe3LSkpIYQQ8ujRIzJixAgiEonIxIkTSV5ennk6biDaxhkREUGmT59OCCHk5cuX7GcMd+7csZqHASG2O18JoXOWzlkFdM5WXqgQVwYyMjJIv379yJAhQ0hSUhIhRHkxLVu2jISEhLArPubGsXnzZvLgwQPzdNgIMJNow4YNpGXLlmT58uWC7W7fvk32799P7t27V5Hd0xvVm5sqhvy2hBDy8ccfk2rVqpGbN2+artNloLzjDA4OVjNJ6DqmuTDmb2oN81Xf34HOWTpn6Zyt/NDABhUSExMRGxuLpKQktc+YPDxVqlTBG2+8gfv377OlPpjwbi8vL1SpUoV1sLSzswMADB8+HLVq1aqIIeiNPmNlkMlkAICBAwciOjoa+/fvR1xcHADg6tWrABRpDOrVq4fevXujTp06Ju694aSlpfHSKnBD0Q39bZl9ly5disTERMTExFTUMHRijHH6+PioOQlboqO7MX9TwLLnK6DfeBkqw5xNTU1Fbm4u+76yzlljjNNa5qwxf1PA8uesqaFCXClSqRTjx49HkyZNMHr0aDRs2BDnzp0DoLzI7O3tUVRUhB07dmD06NFo1KgRdu7ciRMnTrDHef78Ofz9/REeHm6WceiDvmOVSqX49ddf2fdyuRyenp4YNGgQ5HI5FixYgM6dO6NZs2bIzMy02NBuqVSKcePGoW3btnjjjTcwatQotf4a+tuqRkdZAqYYp6ViS2MF9B9vZZmzJSUlGDNmDFq0aIEuXbpg2LBhSE9Pr3Rz1hTjtFRsaawVirlVgZZAbm4uefPNN0mnTp3I1atXyf3790m3bt1Ihw4deO1WrlxJfHx8SN++fQkhhNy4cYMMGzaMODo6kokTJ5Jx48YRDw8PsmbNGkKIZaqyDR3rW2+9xfrSMDx79ozUqFGDiEQiMnToUJKSklKBIzCMjIwM0qVLF9KpUydy9uxZ8tNPP5HGjRuTNm3akPv377PtrP23tZVxEmJbYyXE8PFa+5yVSqVk2LBhpFWrVuTkyZNk+fLlpH79+qRdu3bk7t27bDtr/31tZZyE2NZYKxoqxBFCLl68SGrWrEn+/fdfdtv69evJm2++yV4kP/74I4mIiCDbtm3j+R3I5XKyePFiMnbsWNKrVy9y7ty5Cu+/IRg6VtVJcvz4ceLu7k4aNWpErly5UqF9LwuHDx8m9evX5z3s7t69S8RiMZk6dSrJzMwkGzduJGFhYVb929rKOAmxrbESYvh4rX3OJiQkkJo1a5ItW7aw25KTk0nVqlXJlClTSEZGRqX4fW1lnITY1lgrGirEEULOnDlDRCIRe3GkpaWRRo0akQkTJpC1a9cSQhSh6fn5+bz9rHEVUNaxMrx69Yps3769wvpbXn799Vfi7e3N23bu3Dni4+NDatasSQ4cOEDkcrlahJq1/ba2Mk5CbGushJR9vAzWNmevXbtGXFxcSFxcHCGEsFGVP/74I6lZsyb5+++/iVwut/r7sa2MkxDbGmtFYznOARXE4sWLMW/ePOzYsYPd1q5dO3Ts2BGjRo1Cz549ERgYiKCgIDg6OuKzzz7DoEGDcPv2bbi6uvKyXVui0ygXY44VUJTh8fX1ZbNmWxpC4w0LC4OPjw+++uordtvPP/+MMWPGQC6XY9++fRCJRHBxceEdy5J/W1sZJ2BbYwWMO17A8ufswYMHAfBLfNWuXRvBwcHYunUrAKUP2+TJk+Hl5YXdu3dDIpHA1dWVdyxL/n1tZZyAbY3VIjCjAFmhXLx4kYSFhZEmTZqQnj17Eg8PD/LWW2+xJorc3FwSFxdH2rRpw0v8eP36dVK9enU2WaY1YEtjJUR4vP379yeJiYmkqKiIfPXVV0QkEpE2bdoQd3d3Ur9+fSKVSskPP/xAqlatau7u642tjJMQ2xorIbY33v3795OqVavyrAKMCa2goIB88sknpGbNmuTly5eEEEIKCwsJIYRs2bKFeHl5se8tHVsZJyG2NVZLwmaEuBkzZpDevXsTQhQX1s2bN0l4eDiZOHEiSU5OJoQQcvnyZVK7dm2SmprKqnGlUilbvsRasKWxEqJ5vBMmTCCpqamEEEL+/fdf8sMPP/DKuCxdupS0a9eOZGVlmaXfhmIr4yTEtsZKiG2N98yZM6RHjx7kgw8+ID179iTNmjVTa3Ps2DHSvHlzMmnSJEKI0qx24sQJEhAQQG7cuFGhfS4LtjJOQmxrrJZGpRfi5HI5ycrKIu3atSMzZ84khChXB6tXryZNmzYlK1asIIQQcv/+fSISiUhsbCy7/59//kmaNGlCrl69WvGdNxBbGishho1XFYlEQvr160emTJlSYf0tK7YyTkJsa6yE2NZ4mYf2w4cPyfLly8mTJ0/IlStXiKurK/n5558JIcqEroWFheS7774jbm5uZM+ePWyViUWLFpGOHTtatK+UrYyTENsaq6VSKYW42NhYtZVps2bNyPjx4wkhSqfK4uJiMmDAANKvXz/y7Nkzkp+fT4YMGUJcXV3JhAkTyIgRI4iHhweZO3euxV5gtjRWQgwfb//+/cmTJ0/Ytvfv3ycPHz4kI0aMIJGRkeTChQsV13kDsJVxEmJbYyWEjpcQZZZ9qVRKPvroI+Lv78+Om/ksJyeHfPLJJ8TDw4N06NCBDBo0iLi4uLCF3C3tPmUr4yTEtsZq6VQqIe6PP/4g1apVIzVq1CBhYWFk7ty55Pnz54QQRf4Zd3d3NvqFWQXs3r2bVKtWjS0YnZ+fTz755BMycuRIMmLECIst42FLYyWk7OMNDQ3lhaQvW7aM1KhRg7Rv316tiLIlYCvjJMS2xkoIHe/cuXNZdw65XM4+sJ88eUJCQ0PZepiqpaN27dpF5s2bRyZMmGCRpcFsZZyE2NZYrYVKI8RdvnyZ1KlTh6xYsYLcuHGDrF69mvj7+5OJEyeSrKwsNtkls9rlFsn19fVlVb8MjArYErGlsRJS/vFu2LCBfZ+cnMwzIVsStjJOQmxrrITQ8XLHm56eTghRamfkcjlZvXo1sbe3ZzWOEomEZGdnm63/+mIr4yTEtsZqTVi9EMdI/mvWrCHVqlXjXSQ//vgjadGiBVmyZAkhhJBVq1YROzs7curUKbbN48ePSY0aNcju3bsrtuNlwJbGSojtjNdWxkmIbY2VEDpe1fG2atWKLFy4UG2/9PR00qZNG9K3b18SGxtLunXrRrZs2WKx5jVbGSchtjVWa8TqhTiGTz75hLz++uu8ZIF5eXlk8uTJpFWrVuTBgwdELpeTYcOGkaCgILJgwQJy7do1Mn78eBITE0NevHhhxt4bhi2NlRDbGa+tjJMQ2xorIXS8hCjH26ZNG3L79m1CiFJzQwghGzduJCKRiIjFYtKnTx9SUFBQ4f02FFsZJyG2NVZrwuqEuCNHjpApU6aQFStWkIsXL7Lb9+3bR5ydncnjx48JIcoL6ciRI6RNmzZk+fLlbNspU6aQRo0akaioKNKkSRNy8+bNih2EntjSWAmxnfHayjgJsa2xEkLHy6BtvG3btuWNVyKRkFWrVhGxWEw6dOjACgOWhK2MkxDbGmtlwGqEuKSkJNKnTx8SEBBAhg0bRmJiYoiXlxd7kRUWFpI6deqQcePGEUL4jpSvvfYamThxIvueKSvFrUVoSdjSWAmxnfHayjgJsa2xEkLHa+h4mVxhhBCSkpJCpk2bRn799deKHYQe2Mo4CbGtsVYmrEKIy8/PJ++99x4ZMmQIL9S+efPmZOTIkYQQxapg8+bNRCwWqxXIHTZsGOnUqRP73pJt8rY0VkJsZ7y2Mk5CbGushNDxMpR1vJaKrYyTENsaa2XDKmqnurq6wsnJCSNHjkRkZCRKSkoAAH369MG9e/cAAHZ2dhg8eDD69u2L999/H6dOnQIhBCkpKYiLi8OwYcPY41lyPTZbGitgO+O1lXECtjVWgI63vOO1VGxlnIBtjbXSYTbx0UC4IffMSvXdd98lY8eO5W0rLCwkHTt2JAEBAaRbt24kJCSEtGrViiQkJFR8p8uILY2VENsZr62MkxDbGishdLyEVM7x2so4CbGtsVYmRIQQYm5Bsqy0b98eo0ePxsiRI0EIgVwuh52dHV6+fImbN2/i8uXLiIiIwDvvvGPurpYbWxorYDvjtZVxArY1VoCOt7KO11bGCdjWWK0WMwmP5ebx48ckMDCQXLlyhd3GZDmvbNjSWAmxnfHayjgJsa2xEkLHS0jlHK+tjJMQ2xqrNWMVPnFcSKni8OzZs3B3d0fTpk0BAAsWLMC0adOQmppqzu4ZFVsaK2A747WVcQK2NVaAjreyjtdWxgnY1lgrA/bm7oChME6/ly5dwltvvYWjR49i3LhxKCgowJYtWxAQEGDmHhoPWxorYDvjtZVxArY1VoCOt7KO11bGCdjWWCsFZtMBloPCwkISFRVFRCIRcXJyIkuXLjV3l0yGLY2VENsZr62MkxDbGishdLyVdby2Mk5CbGus1o7VBjZ07doVNWvWxPLly+Hs7Gzu7pgUWxorYDvjtZVxArY1VoCOt7JiK+MEbGus1ozVCnEymQx2dnbm7kaFYEtjBWxnvLYyTsC2xgrQ8VZWbGWcgG2N1ZqxWiGOQqFQKBQKxZaxuuhUCoVCoVAoFAoV4igUCoVCoVCsEirEUSgUCoVCoVghVIijUCgUCoVCsUKoEEehUCgUCoVihVAhjkKhUCgUCsUKoUIchUKxGTp27Ijp06ebuxtqWGq/KBSKZUOFOAqFQhHg5MmTEIlEyMrKMvkx9+zZg4ULFxrtPBQKxTawN3cHKBQKxdbx8fExdxcoFIoVQjVxFAqlUpKfn48RI0bA3d0dwcHBWLZsGe/zrVu3olmzZvDw8EBQUBDeeecdpKamAgDi4+PRqVMnAECVKlUgEokwcuRIAAAhBF9//TWqV68OFxcXNGzYEH/88YfO/mg7pqo5NSIiAosWLWL7Hx4ejn379iEtLQ19+/aFu7s7YmJicOXKFd45zp8/j/bt28PFxQWhoaGYOnUq8vPzy/L1USgUK4AKcRQKpVLy8ccf48SJE/jzzz9x5MgRnDx5ErGxseznxcXFWLhwIW7cuIG9e/fi6dOnrFAVGhqK3bt3AwAePHiA5ORkrFy5EgDw2WefYePGjVizZg3u3LmDDz/8EO+++y5OnTqltT/ajinEd999h7Zt2+LatWvo3bs3hg8fjhEjRuDdd9/F1atXERUVhREjRoCpnHjr1i10794dAwYMwM2bN7Fz506cPXsWH3zwQZm/QwqFYuEQCoVCqWTk5uYSR0dHsmPHDnZbeno6cXFxIdOmTRPc59KlSwQAyc3NJYQQcuLECQKAZGZmsm3y8vKIs7MzOX/+PG/fMWPGkLfffltnv4SOSQghHTp04PUrPDycvPvuu+z75ORkAoB8/vnn7LYLFy4QACQ5OZkQQsjw4cPJuHHjeMc9c+YMEYvFpLCwUGffKBSK9UF94igUSqXj8ePHKC4uRuvWrdltPj4+qF27Nvv+2rVrmD9/Pq5fv46MjAzI5XIAQEJCAqKjowWPe/fuXRQVFaFr16687cXFxWjcuLFRx9CgQQP2dWBgIAAgJiZGbVtqaiqCgoIQGxuLR48eYdu2bWwbQgjkcjmePn2KunXrGrV/FArF/FAhjkKhVDpIqYlRE/n5+ejWrRu6deuGrVu3wt/fHwkJCejevTuKi4s17scIegcOHEDVqlV5nzk5OZW/4xwcHBzY1yKRSOM2pk9yuRzjx4/H1KlT1Y4VFhZm1L5RKBTLgApxFAql0hEVFQUHBwf8999/rACTmZmJhw8fokOHDrh//z5evXqFpUuXIjQ0FADUggQcHR0BADKZjN0WHR0NJycnJCQkoEOHDgb3S+iYxqJJkya4c+cOoqKijH5sCoVimdDABgqFUulwd3fHmDFj8PHHH+P48eO4ffs2Ro4cCbFYccsLCwuDo6MjfvjhBzx58gR//fWXWp628PBwiEQi7N+/H2lpacjLy4OHhwdmzpyJDz/8EL/++iseP36Ma9euYdWqVfj111919kvomMZi1qxZuHDhAiZPnozr168jLi4Of/31F6ZMmWK0c1AoFMuCCnEUCqVS8s0336B9+/Z488030aVLF7Rr1w5NmzYFAPj7+2PTpk34/fffER0djaVLl+Lbb7/l7V+1alUsWLAAs2fPRmBgIBvluXDhQsydOxdLlixB3bp10b17d/z999+IjIzU2SdNxzQGDRo0wKlTpxAXF4fXXnsNjRs3xueff47g4GCjnYNCoVgWIqLLeYRCoVAoFAqFYnFQTRyFQqFQKBSKFUKFOAqFQjESEyZMgLu7u+DfhAkTzN09CoVSyaDmVAqFQjESqampyMnJEfzM09MTAQEBFdwjCoVSmaFCHIVCoVAoFIoVQs2pFAqFQqFQKFYIFeIoFAqFQqFQrBAqxFEoFAqFQqFYIVSIo1AoFAqFQrFCqBBHoVAoFAqFYoVQIY5CoVAoFArFCqFCHIVCoVAoFIoV8n/txCzu50jwjwAAAABJRU5ErkJggg==" + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "from pytesmo.validation_framework.adapters import AnomalyClimAdapter\n", + "ds_mask = AnomalyClimAdapter(ismn_reader, columns=['soil_moisture'], \n", + " timespan=(datetime(1991,1,1), datetime(2020,12,31)))\n", + "anom = ds_mask.read(ids[0])\n", + "pprint(anom)\n", + "anom.plot(title='Anomaly (wrt. 1991-2020 avg.)', ylabel='SM m3m-3', figsize=(7,2))" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "pycharm": { + "name": "#%% md\n" + } + }, + "source": [ + "## ColumnCombineAdapter\n", + "\n", + "This adapter is used to combine multiple columns of a dataset after reading the time series. It takes any function that can be applied accross multiple columns and will add a new column of the chosen name to the data frame.\n", + "E.g. in the following example we combine the `proc_flag` and `snow_prob` and `frozen_prob` column into a new column `good` that is 'True' when all the chosen columns are 0 and otherwise 'False'. Then we apply a `SelfMaskingAdapter` that uses the newly added column to filter the dataset. This is just for demonstration, we could also just use the `AdvancedMaskingAdapter` for this example. A more common use case for the `ColumnCombineAdapter` would be to compute the average when multiple soil moisture fields are available in a dataset.\n", + "\n", + "This example also shows that you can stack multiple adapters together. If they depend on each other, it is important to notice that the innermost adapter will be called first!" ] + }, + { + "cell_type": "code", + "execution_count": 20, + "metadata": { + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " sm proc_flag snow_prob frozen_prob good\n", + "2007-06-03 03:33:38.999992 84.0 0 0 0 True\n", + "2007-06-03 14:51:50.999997 82.0 0 0 0 True\n", + "2007-06-06 02:31:33.000005 76.0 0 0 0 True\n", + "2007-06-06 13:49:46.999999 75.0 0 0 0 True\n", + "2007-06-08 03:30:13.000023 83.0 0 0 0 True\n", + "... ... ... ... ... ...\n", + "2014-06-25 13:51:37.000006 83.0 0 0 0 True\n", + "2014-06-27 03:32:04.999979 79.0 0 0 0 True\n", + "2014-06-27 14:50:18.000019 81.0 0 0 0 True\n", + "2014-06-30 02:30:05.000000 77.0 0 0 0 True\n", + "2014-06-30 13:48:17.999999 83.0 0 0 0 True\n", + "\n", + "[713 rows x 5 columns]\n" + ] + } + ], + "source": [ + "from pytesmo.validation_framework.adapters import ColumnCombineAdapter, SelfMaskingAdapter\n", + "\n", + "def select_good(x):\n", + " return (x['proc_flag'] == 0) and (x['snow_prob'] == 0) and (x['frozen_prob'] == 0)\n", + "\n", + "ds_mask = SelfMaskingAdapter(ColumnCombineAdapter(ascat_reader, func=select_good, new_name='good'), \n", + " '==', True, 'good')\n", + "pprint(ds_mask.read(1814367)[['sm', 'proc_flag', 'snow_prob', 'frozen_prob', 'good']])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { @@ -624,7 +1026,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.9" + "version": "3.9.16" } }, "nbformat": 4, diff --git a/src/pytesmo/time_series/grouping.py b/src/pytesmo/time_series/grouping.py index 1e905f6c..180a8919 100644 --- a/src/pytesmo/time_series/grouping.py +++ b/src/pytesmo/time_series/grouping.py @@ -26,9 +26,6 @@ # NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, # EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. -# Author: Christoph Paulik christoph.paulik@geo.tuwien.ac.at -# Creation date: 2014-06-30 - """ Module provides grouping functions that can be used together with pandas @@ -36,12 +33,16 @@ there are three products per month with timestamps on the 10th 20th and last of the month """ +from dataclasses import dataclass +from typing import Optional, Union, Tuple, List import pandas as pd import numpy as np -from datetime import date +from datetime import date, datetime import calendar +from cadati.conv_doy import doy + def group_by_day_bin(df, bins=[1, 11, 21, 32], start=False, dtindex=None): @@ -153,3 +154,239 @@ def grouped_dates_between(start_date, end_date, bins=[1, 11, 21, 32], start=Fals tstamps = grp.sum().index.to_pydatetime().tolist() return tstamps + + +@dataclass +class YearlessDatetime: + """ + Container class to store Datetime information without a year. This is + used to group data when the year is not relevant (e.g. seasonal analysis). + Only down to second. Used by + :class:`pytesmo.validation_framework.metric_calculators_adapters.TsDistributor` + """ + month: int + + day: int = 1 + hour: int = 0 + minute: int = 0 + second: int = 0 + + @property + def __ly(self): + return 2400 # arbitrary leap year + + def __ge__(self, other: 'YearlessDatetime'): + return self.to_datetime(self.__ly) >= other.to_datetime(self.__ly) + + def __le__(self, other: 'YearlessDatetime'): + return self.to_datetime(self.__ly) <= other.to_datetime(self.__ly) + + def __lt__(self, other: 'YearlessDatetime'): + return self.to_datetime(self.__ly) < other.to_datetime(self.__ly) + + def __gt__(self, other: 'YearlessDatetime'): + return self.to_datetime(self.__ly) > other.to_datetime(self.__ly) + + def __repr__(self): + return f"****-{self.month:02}-{self.day:02}" \ + f"T{self.hour:02}:{self.minute:02}:{self.second:02}" + + @property + def doy(self) -> int: + """ + Get day of year for this date. Assume leap year! + i.e.: 1=Jan.1st, 366=Dec.31st, 60=Feb.29th. + """ + return doy(self.month, self.day, year=None) + + @classmethod + def from_datetime(cls, dt: datetime): + """ + Omit year from passed datetime to create generic datetime. + """ + return cls(dt.month, dt.day, dt.hour, dt.minute, dt.second) + + def to_datetime(self, years: Optional[Union[Tuple[int, ...], int]]) \ + -> Union[datetime, List, None]: + """ + Convert generic datetime to datetime with year. + Feb 29th for non-leap-years will return None + """ + dt = [] + + for year in np.atleast_1d(years): + if not calendar.isleap(year) and self.doy == 60.: + continue + else: + d = datetime(year, self.month, self.day, self.hour, + self.minute, self.second) + dt.append(d) + + if len(dt) == 1: + return dt[0] + elif len(dt) == 0: + return None + else: + return dt + + +class TsDistributor: + + def __init__(self, + dates=None, + date_ranges=None, + yearless_dates=None, + yearless_date_ranges=None): + """ + Build a data distibutor from individual dates, date ranges, generic + dates (without specific year) and generic date ranges. + + Components: + - individual datetime objects for distinct dates + - generic datetime objects for dates without specific a year + - date range / datetime tuple + i.e. ALL datetimes between the 2 passed dates (start, end) + the start date must be earlier than the end date + - generic date range / generic datetime tuple + i.e. ALL datetimes between 2 generic dates (for any year) + + Parameters + ---------- + dates : Tuple[datetime, ...] or Tuple[str, ...] or pd.DatetimeIndex + Individual dates (that also have a year assigned). + date_ranges: Tuple[Tuple[datetime, datetime], ...] + A list of date ranges, consisting of a start and end date for each + range. The start date must be earlier in time than the end date. + yearless_dates: Tuple[YearlessDatetime,...] or Tuple[datetime...] + A list of generic dates (that apply to any year). + Can be passed as a list of + - YearlessDatetime objects + e.g. YearlessDatetime(5,31,12,1,10), ie. May 31st 12:01:10 + - pydatetime objects (years will be ignored, duplicates dropped) + yearless_date_ranges: [Tuple[YearlessDatetime, YearlessDatetime], ...] + A list of generic date ranges (that apply to any year). + """ + + self.dates = dates + self.date_ranges = date_ranges + self.yearless_dates = yearless_dates + self.yearless_date_ranges = yearless_date_ranges + + def __repr__(self): + s = [] + for var in ['dates', 'date_ranges', 'yearless_dates', + 'yearless_date_ranges']: + val = getattr(self, var) + s.append(f"#{var}={len(val) if val is not None else 0}") + + return f"{self.__class__.__name__}({', '.join(s)})" + + def select(self, + df: Union[pd.DataFrame, pd.Series, pd.DatetimeIndex], + set_nan=False): + """ + Select rows from data frame or series with mathing date time indices. + + Parameters + ---------- + df: pd.DataFrame or pd.Series + Must have a date time index, which is then filtered based on the + dates. + set_nan: bool, optional (default: False) + Instead of dropping rows that are not selected, set their values to + nan. + + + Returns + ------- + df: pd.DataFrame or pd.Series + The filterd input data + + """ + if isinstance(df, pd.DatetimeIndex): + idx = df + else: + idx = df.index + + if not isinstance(idx, pd.DatetimeIndex): + raise ValueError(f"Expected a DatetimeIndex, " + f"but got {type(df.index)}.") + + mask = self.filter(idx) + + if set_nan: + df[~mask] = np.nan + return df + else: + return df[mask] + + def filter(self, idx: pd.DatetimeIndex): + """ + Filter datetime index for a TimeSeriesDistributionSet + + Parameters + ---------- + idx: pd.DatetimeIndex + Datetime index to split using the set + + Returns + ------- + idx_filtered: pd.DatetimeIndex + Filtered Index that contains dates for the set + """ + + mask = pd.DataFrame(index=idx.copy()) + + if self.dates is not None: + _idx_dates = idx.intersection(pd.DatetimeIndex(self.dates)) + mask['dates'] = False + mask.loc[_idx_dates, 'dates'] = True + + if self.date_ranges is not None: + for i, drange in enumerate(self.date_ranges): + start, end = drange[0], drange[1] + if start > end: + start, end = end, start + mask[f"range{i}"] = (idx >= start) & (idx <= end) + + if self.yearless_dates is not None: + arrs = np.array([]) + for d in self.yearless_dates: + dts = d.to_datetime(np.unique(idx.year)) + if dts is None: + continue + else: + arrs = np.append(arrs, dts) + _idx_dates = idx.intersection(pd.DatetimeIndex(arrs)) + mask['gen_dates'] = False + mask.loc[_idx_dates, 'gen_dates'] = True + + # avoid loop like: + # cond = ["__index_month == {}".format(m) for m in months] + # selection = dat.query(" | ".join(cond)).index + + if self.yearless_date_ranges is not None: + for i, gdrange in enumerate(self.yearless_date_ranges): + for y in np.unique(idx.year): + + if not calendar.isleap(y) and (gdrange[0].doy == 60): + start = YearlessDatetime(3, 1) + else: + start = gdrange[0] + + if (not calendar.isleap(y)) and (gdrange[1].doy == 60): + end = YearlessDatetime(2, 28, 23, 59, 59) + else: + end = gdrange[1] + + start_dt = start.to_datetime(years=y) + + if end < start: + end_dt = end.to_datetime(years=y + 1) + else: + end_dt = end.to_datetime(years=y) + + mask[f"gen_range{y}-{i}"] = (idx >= start_dt) & ( + idx <= end_dt) + + return mask.any(axis=1, bool_only=True) diff --git a/src/pytesmo/validation_framework/adapters.py b/src/pytesmo/validation_framework/adapters.py index 8229adf5..115d71ba 100644 --- a/src/pytesmo/validation_framework/adapters.py +++ b/src/pytesmo/validation_framework/adapters.py @@ -31,10 +31,8 @@ import operator -import pandas as pd from pytesmo.time_series.anomaly import calc_anomaly from pytesmo.time_series.anomaly import calc_climatology -from pytesmo.utils import deprecated from pandas import DataFrame import numpy as np import warnings @@ -90,9 +88,10 @@ def __init__(self, cls, data_property_name="data", read_name=None): setattr(self, read_name, self._adapt_custom) def __get_dataframe(self, data): - if ((not isinstance(data, DataFrame)) and - (hasattr(data, self.data_property_name)) and - (isinstance(getattr(data, self.data_property_name), DataFrame))): + if (not isinstance(data, DataFrame)) and \ + (hasattr(data, self.data_property_name)) and \ + (isinstance(getattr(data, self.data_property_name), + DataFrame)): data = getattr(data, self.data_property_name) return data @@ -131,14 +130,14 @@ def grid(self): return self.cls.grid -@deprecated("`MaskingAdapter` is deprecated, use `SelfMaskingAdapter` " - "or `AdvancedMaskingAdapter` instead.") class MaskingAdapter(BasicAdapter): """ Transform the given class to return a boolean dataset given the operator and threshold. This class calls the read_ts and read methods of the given instance and applies boolean masking to the returned data - using the given operator and threshold. + using the given operator and threshold. This adapter does not filter the + time series (see the AdvancedMaskingAdapter and SelfMaskingAdapter for + that) but only turns it into a boolean dataset. Parameters ---------- @@ -267,8 +266,8 @@ class AdvancedMaskingAdapter(BasicAdapter): Reader object, has to have a `read_ts` or `read` method or a method name must be specified in the `read_name` kwarg. The same method will be available for the adapted version of the reader. - filter_list: list[tuple, tuple, tuple] - [column_name, operator, threshold] + filter_list: list[tuple] + [(column_name, operator, threshold), ...] 'column_name': string name of the column to apply the operator to 'operator': Callable or str; @@ -555,16 +554,17 @@ class TimestampAdapter(BasicAdapter): name must be specified in the `read_name` kwarg. The same method will be available for the adapted version of the reader. time_offset_fields: str, list or None - name or list of names of the fields that provide information on the time offset. - If a list is given, all values will contribute to the offset, assuming that - each refers to the previous. For instance: + name or list of names of the fields that provide information on the + time offset. + If a list is given, all values will contribute to the offset, assuming + that each refers to the previous. For instance: offset = minutes + seconds in the minute + µs in the second NOTE: np.nan values are counted as 0 offset NOTE: if None, no offset is considered time_units: str or list - time units that the time_offset_fields are specified in. If a list is given, - it should have the same size as the 'time_offset_fields' parameter - Can be any of the np.datetime[64] units: + time units that the time_offset_fields are specified in. If a list is + given, it should have the same size as the 'time_offset_fields' + parameter. Can be any of the np.datetime[64] units: https://numpy.org/doc/stable/reference/arrays.datetime.html base_time_field: str, optional. Default is None. If a name is provided, the generic time field will be searched for @@ -573,19 +573,20 @@ class TimestampAdapter(BasicAdapter): base_time_reference: str, optional. Default is None. String of format 'YYYY-mm-dd' that can be specified to tranform the 'base_time_field' from [units since base_time_reference] to - np.datetime[64]. If not provided, it will be assumed that the base_time_field - is already in np.datetime[64] units + np.datetime[64]. If not provided, it will be assumed that the + base_time_field is already in np.datetime[64] units base_time_units: str, optional. Default is "D" - Units that the base_time_field is specified in. Only applicable with 'base_time_reference' + Units that the base_time_field is specified in. Only applicable with + 'base_time_reference' replace_index: bool, optional. Default is True. If True, the exact timestamp is used as index. Else, it will be added to the dataframe on the column 'output_field' output_field: str, optional. Default is None. - If a name is specified, an additional column is generated under the name, - with the exact timestamp. Only with 'replace_index' == False + If a name is specified, an additional column is generated under the + name, with the exact timestamp. Only with 'replace_index' == False drop_original: bool, optional. Default is True. - Whether the base_time_field and time_offset_fields should be dropped in the - final DataFrame + Whether the base_time_field and time_offset_fields should be dropped + in the final DataFrame """ def __init__(self, @@ -617,12 +618,13 @@ def __init__(self, if not replace_index and output_field is None: raise ValueError( "'output_field' should be specified in case the new timestamp" - "should not be used as index. Alternatively, set 'replace_index' to True" + "should not be used as index. Alternatively, set " + "'replace_index' to True" ) elif replace_index and output_field is not None: warnings.warn( - "Ignoring the 'output_field' value. Set 'replace_index' to True to" - "avoid this behavior") + "Ignoring the 'output_field' value. Set 'replace_index' to " + "True to avoid this behavior") else: self.output_field = output_field @@ -639,7 +641,9 @@ def convert_generic(self, return time_date def add_offset_cumulative(self, data: DataFrame) -> np.array: - """Return an array of timedelta calculated with all the time_offset_fields""" + """ + Return an array of timedelta calculated with all the time_offset_fields + """ total_offset = np.full(data.index.shape, 0, dtype='timedelta64[s]') for field, unit in zip(self.time_offset_fields, self.time_units): total_offset += data[field].map( diff --git a/src/pytesmo/validation_framework/metric_calculators_adapters.py b/src/pytesmo/validation_framework/metric_calculators_adapters.py index dd29a9e2..d6dad14a 100644 --- a/src/pytesmo/validation_framework/metric_calculators_adapters.py +++ b/src/pytesmo/validation_framework/metric_calculators_adapters.py @@ -24,25 +24,31 @@ # ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT # (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF # THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. - - """ Metric Calculator Adapters change how metrics are calculated by calling the `calc_metric` function of the adapted calculator instead of the unadapted version. """ +from pytesmo.time_series.grouping import YearlessDatetime, TsDistributor from pytesmo.validation_framework.metric_calculators import ( - PairwiseIntercomparisonMetrics, - TripleCollocationMetrics -) + PairwiseIntercomparisonMetrics, TripleCollocationMetrics) import warnings import numpy as np +from cadati.conv_doy import days_past -class MonthsMetricsAdapter(object): +def days_in_month(month: int) -> int: """ - Adapt MetricCalculators to calculate metrics for groups across months + Get number of days in this month (in a LEAP YEAR) + """ + return days_past[month] - days_past[month - 1] + + +class SubsetsMetricsAdapter: + """ + Adapt MetricCalculators to calculate metrics for groups of temporal + subsets (also across multiple years). """ _supported_metric_calculators = ( @@ -50,101 +56,64 @@ class MonthsMetricsAdapter(object): TripleCollocationMetrics, ) - def __init__(self, calculator, sets=None): + def __init__(self, calculator, subsets, group_results='tuple'): """ Add functionality to a metric calculator to calculate validation - metrics for subsets of certain months in a time series (e.g. seasonal). + metrics for subsets of certain datetimes in a time series + (e.g. seasonal). Parameters ---------- calculator : PairwiseIntercomparisonMetrics or TripleCollocationMetrics A metric calculator to adapt. Preferably an instance of a metric calculator listed in `_supported_metric_calculators` - sets : dict, optional (default: None) - Define groups of data. With group names as key and a list of - months (1-12) that belong to the group as values. - - e.g. {'Group1': [4,5,6,7,8,9], 'Group2': [10,11,12,1,2,3]} will - split the data used by the metric calculator into 2 groups. - One using only observations made between April and September, - and one using observations from the rest of the year. - - The name will be used in the results to distinguish between the - same metrics for different groups: - e.g. ('Group1', 'BIAS'): ..., ('Group2', 'BIAS'): ..., etc. - - The default groups are based on 4 seasons plus one group that uses - all data (as the unadapted metric calculator would do): - {'DJF': [12,1,2], 'MAM': [3,4,5], 'JJA': [6, 7, 8], - 'SON': [9, 10, 11], 'ALL': list(range(1, 13))} + subsets : dict[str, TsDistributor], optional (default: None) + Define subsets of data. With group names as key and a + data distributor as values. + group_results: str, optional (default: 'tuple') + How to group the results. + - 'tuple' will group the results by (group, metric) + - 'join' will join group and metric name with a '|' """ if not isinstance(calculator, self._supported_metric_calculators): warnings.warn(f"Adapting {calculator.__class__} is not supported.") + self.cls = calculator - if sets is None: - sets = { - "DJF": [12, 1, 2], - "MAM": [3, 4, 5], - "JJA": [6, 7, 8], - "SON": [9, 10, 11], - "ALL": list(range(1, 13)), - } + self.subsets = subsets + self.group_results = group_results - self.sets = sets + assert group_results in ('tuple', 'join'), \ + f"Unknown group_results: {group_results}" # metadata metrics and lon, lat, gpi are excluded from applying # seasonally self.non_seas_metrics = ["gpi", "lon", "lat"] - if self.cls.metadata_template is not None: - self.non_seas_metrics += list(self.cls.metadata_template.keys()) + if hasattr(self.cls, 'metadata_template'): + if self.cls.metadata_template is not None: + self.non_seas_metrics += list( + self.cls.metadata_template.keys()) all_metrics = calculator.result_template subset_metrics = {} # for each subset create a copy of the metric template - for name in sets.keys(): + for name in subsets.keys(): for k, v in all_metrics.items(): - if k in self.non_seas_metrics: - subset_metrics[k] = np.array(v) - else: - subset_metrics[(name, k)] = np.array(v) + subset_metrics[self._genname(name, k)] = np.array(v) self.result_template = subset_metrics - @staticmethod - def filter_months(df, months, dropna=False): - """ - Select only entries of a time series that are within certain month(s) - - Parameters - ---------- - df : pd.DataFrame - Time series (index.month must exist) that is split up into the - selected groups. - months : list - Months for which data is kept, e.g. [12,1,2] to keep data for - winter - dropna : bool, optional (default: False) - Drop lines for months that are not to be kept, if this is false, - the original index is not changed, but filtered values are replaced - with nan. - - Returns - ------- - df_filtered : pd.DataFrame - The filtered series - """ - dat = df.copy(True) - dat["__index_month"] = dat.index.month - cond = ["__index_month == {}".format(m) for m in months] - selection = dat.query(" | ".join(cond)).index - dat.drop("__index_month", axis=1, inplace=True) - - if dropna: - return dat.loc[selection] + def _genname(self, setname: str, metric: str) -> str or tuple: + if metric in self.non_seas_metrics: + k = f"{metric}" + elif self.group_results == 'tuple': + k = (f"{setname}", *np.atleast_1d(metric)) + elif self.group_results == 'join': + k = f"{setname}|{metric}" else: - dat.loc[dat.index.difference(selection)] = np.nan - return dat + raise NotImplementedError( + f"Unknown group_results: {self.group_results}") + return k def calc_metrics(self, data, gpi_info): """ @@ -161,14 +130,63 @@ def calc_metrics(self, data, gpi_info): """ dataset = self.result_template.copy() - for setname, months in self.sets.items(): - df = self.filter_months(data, months=months, dropna=True) + for setname, distr in self.subsets.items(): + df = distr.select(data) ds = self.cls.calc_metrics(df, gpi_info=gpi_info) for metric, res in ds.items(): - if metric in self.non_seas_metrics: - k = f"{metric}" - else: - k = (f"{setname}", *np.atleast_1d(metric)) + k = self._genname(setname, metric) dataset[k] = res return dataset + + +class MonthsMetricsAdapter(SubsetsMetricsAdapter): + """ + Adapt MetricCalculators to calculate metrics for groups across months + """ + + def __init__(self, calculator, sets=None): + """ + Add functionality to a metric calculator to calculate validation + metrics for subsets of certain months in a time series (e.g. seasonal). + + Parameters + ---------- + calculator : PairwiseIntercomparisonMetrics or TripleCollocationMetrics + A metric calculator to adapt. Preferably an instance of a metric + calculator listed in `_supported_metric_calculators` + sets : dict, optional (default: None) + Define groups of data. With group names as key and a list of + months (1-12) that belong to the group as values. + + e.g. {'Group1': [4,5,6,7,8,9], 'Group2': [10,11,12,1,2,3]} will + split the data used by the metric calculator into 2 groups. + One using only observations made between April and September, + and one using observations from the rest of the year. + + The name will be used in the results to distinguish between the + same metrics for different groups: + e.g. ('Group1', 'BIAS'): ..., ('Group2', 'BIAS'): ..., etc. + + The default groups are based on 4 seasons plus one group that uses + all data (as the unadapted metric calculator would do): + {'DJF': [12,1,2], 'MAM': [3,4,5], 'JJA': [6, 7, 8], + 'SON': [9, 10, 11], 'ALL': list(range(1, 13))} + """ + if sets is None: + sets = { + 'DJF': [12, 1, 2], + 'MAM': [3, 4, 5], + 'JJA': [6, 7, 8], + 'SON': [9, 10, 11], + 'ALL': list(range(1, 13)), + } + + for name, months in sets.items(): + distr = TsDistributor(yearless_date_ranges=[( + YearlessDatetime(m, 1, 0, 0, 0), + YearlessDatetime(m, days_in_month(m), 23, 59, 59)) + for m in months]) + sets[name] = distr + + super().__init__(calculator, subsets=sets) diff --git a/src/pytesmo/validation_framework/validation.py b/src/pytesmo/validation_framework/validation.py index 72d37592..cacfe414 100644 --- a/src/pytesmo/validation_framework/validation.py +++ b/src/pytesmo/validation_framework/validation.py @@ -98,7 +98,7 @@ class Validation(object): method of these datasets has to return pandas.DataFrames with only boolean columns. True means that the observations at this timestamp should be masked and False means that it should be kept. - scaling : string, None or class instance + scaling : str or None or class instance - If set then the data will be scaled into the reference space using the method specified by the string using the :py:class:`pytesmo.validation_framework.data_scalers.DefaultScaler` @@ -262,7 +262,8 @@ def calc( ) except Exception as e: raise eh.DataManagerError( - f"Getting the data for gpi {gpi_info} failed!") + f"Getting the data for gpi {gpi_info} failed with" + f" error: {e}") # if no data is available continue with the next gpi if len(df_dict) == 0: @@ -453,7 +454,7 @@ def dummy_result(): except Exception as e: raise eh.ScalingError( f"Scaling failed for {result_key} for gpi" - f" {gpi_info}!" + f" {gpi_info} with error {e}!" ) # Drop the scaling reference if it was not in the intended diff --git a/tests/test_validation_framework/test_metric_calculators_adapters.py b/tests/test_validation_framework/test_metric_calculators_adapters.py new file mode 100644 index 00000000..a442b8fc --- /dev/null +++ b/tests/test_validation_framework/test_metric_calculators_adapters.py @@ -0,0 +1,131 @@ +import unittest + +import pandas as pd +import numpy as np + +from pytesmo.time_series.grouping import YearlessDatetime, TsDistributor +from datetime import datetime + + +class Test_YearlessDateTime(unittest.TestCase): + + def setUp(self) -> None: + self.past = datetime(1900, 1, 2, 3, 4, 5) + self.future = datetime(2104, 6, 7, 8, 9, 10) + self.yearless = YearlessDatetime(3, 10, 0, 0, 0) + + def test_comparisons(self): + assert self.yearless > YearlessDatetime.from_datetime(self.past) + assert self.yearless < YearlessDatetime.from_datetime(self.future) + assert self.yearless == self.yearless + + def test_doy(self): + assert YearlessDatetime.from_datetime(self.future).doy == 159 + assert YearlessDatetime(12, 31).doy == 366 + assert YearlessDatetime(2, 29).doy == 60 + assert YearlessDatetime(1, 1).doy == 1 + + def test_to_dt(self): + assert YearlessDatetime.from_datetime(self.past).to_datetime( + self.past.year) == self.past + assert YearlessDatetime.from_datetime( + self.future).to_datetime(years=[2104, 2111])[0] == self.future + + +class Test_TimeSeriesDistributionSet(unittest.TestCase): + + def setUp(self) -> None: + df = pd.DataFrame( + index=pd.date_range('2000-01-01T12', '2009-12-31T12', freq='D')) + df['data'] = np.random.rand(df.index.size) + df.loc[np.isin(df.index.month, [1, 7])] = np.nan + df = df.dropna() + self.df = df + + def test_filter_dates_only(self): + dates = ( + datetime(2005, 6, 6, 12), + datetime(2005, 5, 5, 12), + datetime(2005, 4, 4, 12), + datetime(2005, 3, 3, 12), + datetime(2005, 2, 2, 12), + datetime(2005, 2, 2, 1), # not in input/output !! + datetime(2005, 1, 5, 12), # not in input/output !! + ) + + set1 = TsDistributor(dates=dates) + + d = set1.select(self.df) + assert len(d.index) == 5 + assert np.all(dt in d.index for dt in dates[:5]) + + def test_filter_daterange_only(self): + set2 = TsDistributor(date_ranges=[ + (datetime(2004, 12, 20), datetime(2005, 2, 10, 11)), + (datetime(2005, 2, 27), datetime(2005, 3, 1, 12)), + ]) + + d = set2.select(self.df) + assert datetime(2005, 2, 1, 12) in d.index + assert len(d.index) == (12 + (10 - 1)) + 3 + + def test_filter_yearless_dates_only(self): + yearless_dates = ( + YearlessDatetime(6, 6, 12, 0, 0), + YearlessDatetime(2, 29, 12, 0, 0), + YearlessDatetime.from_datetime(datetime(2000, 5, 5, 12)), + YearlessDatetime(1, 1, 12), # not in input/output !! + ) + set3 = TsDistributor(yearless_dates=yearless_dates) + + d = set3.select(self.df) + assert datetime(2005, 5, 5, 12) in d.index + assert datetime(2008, 2, 29, 12) in d.index + assert len(d.index) == 2 * len(np.unique(self.df.index.year)) + 3 + + def test_filter_yearless_date_ranges_only(self): + set4 = TsDistributor(yearless_date_ranges=[ + (YearlessDatetime(12, 20), + YearlessDatetime(2, 10, 0)), # 12 + 9 elements + (YearlessDatetime(2, 27), YearlessDatetime(2, 29, 12)) # 3 or 2 + ]) + d = set4.select(self.df) + ny = len(np.unique(self.df.index.year)) + assert datetime(2007, 12, 21, 12) in d.index + assert datetime(2004, 2, 29, 12) in d.index + assert len(d.index) == (12 * ny) + (9 * (ny - 1)) + (3 * ny - 7) + + def test_filter_all_in_one(self): + dates = ( + datetime(2005, 4, 4, 12), + datetime(2005, 5, 5, 12), + ) + date_ranges = [(datetime(2004, 4, 6), datetime(2004, 4, 8, 12))] + yearless_dates = [YearlessDatetime(4, 10, 12)] + yearless_date_ranges = [(YearlessDatetime(2, 27), + YearlessDatetime(2, 29, 23))] + + set = TsDistributor( + dates=dates, + yearless_dates=yearless_dates, + date_ranges=date_ranges, + yearless_date_ranges=yearless_date_ranges, + ) + d = set.select(self.df) + + assert np.all(dt in d.index for dt in dates) + for dt in [ + datetime(2004, 4, 6, 12), + datetime(2004, 4, 7, 12), + datetime(2004, 4, 8, 12) + ]: + assert dt in d.index + + assert datetime(2008, 4, 10, 12) in d.index + + for dt in [ + datetime(2007, 2, 27, 12), + datetime(2007, 2, 28, 12), + datetime(2008, 2, 29, 12) + ]: + assert dt in d.index