{ "cells": [ { "cell_type": "markdown", "id": "0dd05fd2", "metadata": {}, "source": [ "# Creating new *unseen-awg* weather generators" ] }, { "cell_type": "markdown", "id": "925acd71", "metadata": {}, "source": [ "The provided datasets described in {ref}`data-download` can be used to create new *unseen-awg* weather generators. Apart from the reforecast-based weather generators that are our main focus, it is also possible to create weather generators based on the ERA5 reanalysis." ] }, { "cell_type": "markdown", "id": "41abdcd4", "metadata": {}, "source": [ "## Using Snakemake to configure weather generators\n", "\n", "Workflows for creating weather generators from raw reforecast data (or ERA5 data) are specified in the `workflow/` directory. Using the provided data ({ref}`data-download`), it is possible to avoid using raw data and to skip the first preprocessing steps in which data is regridded, aggregated to a daily scale and unified.\n", "\n", "Using the [Snakemake](https://snakemake.github.io/) workflow manager, creating a weather generator is as easy as:\n", "\n", "```bash\n", "conda activate unseen-awg # activate conda environment\n", "cd # navigate to unseen-awg base directory\n", "snakemake --dry-run --cores 1 --configfile configs/weather_generators/reforecasts.yaml\n", "```\n", "\n", "This loads the configuration file `configs/weather_generators/reforecasts.yaml` that specifies the weather generator's parameters and produces:\n", "- a weather generator instance\n", "- a small artificial weather dataset generated with the instance\n", "- a tuning run to find ideal values for some of *unseen-awg*'s parameters\n", "- and if necessary, it also generates a post-processed dataset of impact-relevant variables (e.g., including bias correction). \n", "\n", "Adding the `--dry-run` flag allows testing the command without actually producing output, while `--cores 1` indicates that a single core should be used (which is justified here, given it's a dry run). If other outputs are desired, one can vary the configuration file or adjust the `all` rule defined in `workflow/Snakefile`.\n", "\n", "### Configuration files\n", "\n", "We collected configurations for different weather generator setups in `configs/weather_generators/`. In our study, we use:\n", "- `configs/weather_generators/reforecasts.yaml`\n", "- `configs/weather_generators/reforecasts_no_bias_correction.yaml`\n", "- `configs/weather_generators/era5.yaml`\n", "\n", "These `.yaml` files specify parameter values for all necessary preprocessing steps, some of which are described in more detail in {ref}`preprocessing`.\n", "\n", "Parameters for the rules that define how to preprocess raw data cannot be changed when working with the provided preprocessed datasets ({ref}`data-download`):\n", "- `preprocess_era5_single`\n", "- `preprocess_circulation_era5`\n", "- `preprocess_impact_variables_era5`\n", "- `preprocess_circulation_reforecasts`\n", "- `preprocess_impact_variables_reforecasts`\n", "\n", "If the downloaded data is stored using the correct paths (described in {ref}`data-download`), Snakemake will recognize it as the output of preprocessing rules and skip the early steps of the pipeline.\n", "\n", "### A note on resource demands\n", "\n", "The Snakemake workflow is computationally expensive, especially if the large reforecast dataset is used. While we tried to keep requirements manageable, creating the full *unseen-awg* weather generator will most likely only be possible on HPC systems, given the demands for disk space, memory and compute. Running smaller test cases (e.g., ERA5-based *unseen-awg* weather generators) should also be possible on personal devices.\n", "\n", "Specifying resources (memory, CPU cores, runtime, etc.) on a per-rule basis is possible, and Snakemake also allows creating system-specific profiles. The `resources` section of each Snakemake rule defined in `workflow/rules` will likely have to be adapted to the system the workflow is executed with. We used the `snakemake-executor-plugin-slurm`, which allows Snakemake to execute rules as jobs on a Slurm cluster and created a profile that is adapted to the requirements of our system (not included).\n", "\n", "### Visualizing the Snakemake workflow\n", "The command\n", "```bash\n", "snakemake --configfile configs/weather_generators/reforecasts.yaml --rulegraph | dot -Tpng > rulegraph.png\n", "```\n", "can be used to create a helpful visualization of the rules being executed for the given configuration file (in this case `configs/weather_generators/reforecasts.yaml`) and the current state of the `Snakefile`.\n", "\n", "The resulting rule graph shows how Snakemake produces an *unseen-awg* instance, generates a test simulation with it, performs a tuning run for some of its parameters, and prepares a bias-corrected dataset of impact-relevant variables:\n", "\n", "![Rulegraph of the Snakemake workflow for reforecasts.yaml](images/rulegraph.png)\n", "\n", "The rulegraph includes the \"raw-data preprocessing rules\" that will not be executed if the archived data is used properly. Without bias correction and with a smaller set of targeted outputs, the corresponding graph is much simpler. " ] }, { "cell_type": "markdown", "id": "12324082", "metadata": {}, "source": [ "(preprocessing)=\n", "## Preprocessing details and motivation" ] }, { "cell_type": "markdown", "id": "37a1d7aa", "metadata": {}, "source": [ "Two main lines of work are being executed in the Snakemake pipeline:\n", "- if necessary, impact-relevant variables are bias corrected. To do so, biases are computed, merged, and a Zarr store with corrected data is produced. \n", "- the atmospheric circulation variables are brought into a shape that allows parallelizing similarity computations efficiently.\n", "\n", "`rechunk_ds` is a helper rule that rechunks a Zarr store. For some operations (e.g., computing biases), it is far more efficient to have the entire temporal extent in a single chunk and instead split along a spatial dimension, while for others (e.g., computing similarities) it is more efficient to have the full spatial extent in the same chunk but to split along the time dimension. By default, the archived datasets use the latter chunking option.\n", "\n", "Below, we explain the preprocessing steps necessary for computing similarities efficiently in more detail." ] }, { "cell_type": "code", "execution_count": 15, "id": "4739f5d8", "metadata": { "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "import os\n", "\n", "import matplotlib.colors as mcolors\n", "import matplotlib.pyplot as plt\n", "import numpy as np\n", "import xarray as xr\n", "import yaml\n" ] }, { "cell_type": "code", "execution_count": 16, "id": "cfc8a03f", "metadata": { "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "single_forecast = xr.open_dataset\n", "with open(\"../../configs/paths.yaml\") as file:\n", " paths = yaml.safe_load(file)[\"paths\"]\n", "\n", "path_circulation_ds = os.path.join(\n", " paths[\"dir_preprocessed_datasets\"],\n", " \"preprocessed_circulation_reforecasts/combined_5e06172f.zarr\",\n", ")\n", "\n", "path_restructured = os.path.join(\n", " paths[\"dir_preprocessed_datasets\"],\n", " \"preprocessed_circulation_reforecasts/restructured_5e06172f_f40e9460.nc\",\n", ")\n", "\n", "path_dayofyear_year = os.path.join(\n", " paths[\"dir_preprocessed_datasets\"],\n", " \"preprocessed_circulation_reforecasts/data_year_dayofyear_5e06172f_f40e9460.zarr\",\n", ")\n", "\n", "\n", "path_circulation_ds_era5 = os.path.join(\n", " paths[\"dir_preprocessed_datasets\"],\n", " \"preprocessed_circulation_era5/combined_f3d1f2f7.zarr\",\n", ")\n" ] }, { "cell_type": "markdown", "id": "26d648aa", "metadata": {}, "source": [ "### From initiazlization time and lead time to valid time and sample" ] }, { "cell_type": "markdown", "id": "b91711f6", "metadata": {}, "source": [ "After the initial preprocessing steps, the atmospheric circulation dataset has temporal dimensions `lead_time` and `init_time`, i.e., forecast lead time and initialization time:" ] }, { "cell_type": "code", "execution_count": 17, "id": "27aa7425", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "FrozenMappingWarningOnValuesAccess({'ensemble_member': 11, 'init_time': 2560, 'lead_time': 46, 'latitude': 18, 'longitude': 49})" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds_combined = xr.open_zarr(path_circulation_ds)\n", "ds_combined.dims" ] }, { "cell_type": "markdown", "id": "120ce879", "metadata": {}, "source": [ "The data is spaced unevenly along the `init_time` dimension in intervals of three, four, or less days, because ECMWF ran reforecasts twice weekly in the IFS configuration we downloaded data for.\n", "\n", "We want to compute similarities between atmospheric states with similar dates, and in this case we're interested in the *valid time*, i.e., the actual time the forecast making predictions for. It can be computed using *init_time* + *lead_time*.\n", "\n", "Efficiently parallelizing the similarity computations therefore would require having `valid_time` as a dimension." ] }, { "cell_type": "code", "execution_count": 18, "id": "340d162f", "metadata": { "tags": [ "remove-cell" ] }, "outputs": [], "source": [ "ds_combined = ds_combined.assign_coords(\n", " valid_time=ds_combined.init_time + ds_combined.lead_time\n", ")\n", "unique_valid_times = np.unique(ds_combined.valid_time)\n", "\n", "vt_min = np.datetime64(\"2010-01-01\", \"ns\")\n", "vt_max = np.datetime64(\"2010-02-01\", \"ns\")\n", "\n", "range_vt = np.arange(vt_min, vt_max, np.timedelta64(1, \"D\"))" ] }, { "cell_type": "code", "execution_count": 19, "id": "72c5b738", "metadata": { "tags": [ "remove-cell" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (lead_time: 46, valid_time: 31)> Size: 1kB\n",
       "array([[ True, False, False, ...,  True, False, False],\n",
       "       [False,  True, False, ..., False,  True, False],\n",
       "       [False, False,  True, ..., False, False,  True],\n",
       "       ...,\n",
       "       [False,  True, False, ..., False,  True, False],\n",
       "       [False, False,  True, ..., False, False,  True],\n",
       "       [False, False, False, ..., False, False, False]])\n",
       "Coordinates:\n",
       "  * lead_time   (lead_time) float64 368B 0.0 1.0 2.0 3.0 ... 42.0 43.0 44.0 45.0\n",
       "  * valid_time  (valid_time) datetime64[ns] 248B 2010-01-01 ... 2010-01-31
" ], "text/plain": [ " Size: 1kB\n", "array([[ True, False, False, ..., True, False, False],\n", " [False, True, False, ..., False, True, False],\n", " [False, False, True, ..., False, False, True],\n", " ...,\n", " [False, True, False, ..., False, True, False],\n", " [False, False, True, ..., False, False, True],\n", " [False, False, False, ..., False, False, False]])\n", "Coordinates:\n", " * lead_time (lead_time) float64 368B 0.0 1.0 2.0 3.0 ... 42.0 43.0 44.0 45.0\n", " * valid_time (valid_time) datetime64[ns] 248B 2010-01-01 ... 2010-01-31" ] }, "execution_count": 19, "metadata": {}, "output_type": "execute_result" } ], "source": [ "contains_data = xr.DataArray(\n", " np.full((len(ds_combined.lead_time), len(range_vt)), fill_value=False, dtype=bool),\n", " coords={\"lead_time\": ds_combined.lead_time, \"valid_time\": range_vt},\n", ")\n", "\n", "for lt in ds_combined.lead_time:\n", " for j, vt in enumerate(range_vt):\n", " contains_data.loc[{\"lead_time\": lt, \"valid_time\": vt}] = np.isin(\n", " vt, ds_combined.valid_time.sel(lead_time=lt)\n", " )\n", "contains_data = contains_data.assign_coords(\n", " lead_time=contains_data.lead_time / np.timedelta64(1, \"D\")\n", ")\n", "contains_data.drop_attrs()" ] }, { "cell_type": "markdown", "id": "f90b4543", "metadata": {}, "source": [ "If we would naively use `valid_time` and `lead_time` as coordinates, this would lead to a strongly inflated dataset with a lot of missing data:" ] }, { "cell_type": "code", "execution_count": 20, "id": "cd0062be", "metadata": { "tags": [ "remove-input" ] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm8AAAHWCAYAAAAhG26oAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAARQ5JREFUeJzt3Qd8VFX2wPEDoSSUBFGkg0iR3gQJoNJRRFEQCyiLAq40aSKKSBMp4kpRQAFpumBBbLsKUgRcZePSlCoWEIKAgAIBDC3M/3Oun5n/JLTMJDN5d97v+/k8nXkzeXNf3jBzcu6952bzeDweAQAAgBWyZ3UDAAAAkH4EbwAAABYheAMAALAIwRsAAIBFCN4AAAAsQvAGAABgEYI3AAAAixC8AQAAWCSHRLjz58/Lvn37JH/+/JItW7asbg4AACGhNfePHz8uxYoVk+zZw5ObOXXqlJw5cybTjpcrVy6Jjo7OtONFqogP3jRwK1myZFY3AwCAsEhMTJQSJUqEJXDLXyS/nDt2LtOOGRsbK0WLFjXBZ69evcwGFwZvmnHzvpn1TeF08Qvig/q5hI4JEm42tTXc55gV58n1iIzr4ZZ/I3xmZb6kpCSTrPB+74WaZtw0cLthwg2SPSbjmb7zyedlx4Ad1nxfZ6WID968XaX6RrDhzRAVExXUz2XFudnU1nCfY1acJ9cjMq6HW/6N8JkVOuEeIqSBW0be7wgcExYAAAAsQvAGAABgEYI3AAAAixC8AQAAWITgDQAAwCIEbwAAABaJ+FIhWaXavGpB/dzmzpsl3GxqqxvO0aa2uuEcbWpruM8xK87TpraG+zxTklNC0hY4D5k3AAAAixC8AQAAWITgDQAAwCIEbwAAABYheAMAALAIwRsAAIBFCN4AAAAsQvAGAABgEYr0RkgxSJva6paCqTa11Q3naFNb3XCONrXVls/lpKQkiesRF/Rrwh5k3gAAACxC8AYAAGARgjcAAACLELwBAABYhOANAADAIgRvAAAAFiF4AwAAsAjBGwAAgEVcU6Q3fkG8RMVEBfQzFK4MDZsKCtvUVjeco01tdUsBW5vaGiw3nCPsQuYNAADAIgRvAAAAFiF4AwAAsAjBGwAAgEUI3gAAACxC8AYAAGARgjcAAACLuKbOW0LHBImNjQ3b67mhLpBN52hTW91wjja11Q016Wxqq1ved8G0NSU5JSRtgfOQeQMAALAIwRsAAIBFCN4AAAAsQvAGAABgEYI3AAAAixC8AQAAWITgDQAAwCIEbwAAABZxTJHesWPHyrPPPit9+/aVSZMmmX0ej0dGjhwpM2bMkCNHjki9evVk6tSpUqVKlbC0icKVkVG4UlFINDS4Hs56z3E93F38OCkpSeJ6xAX9mrCHIzJva9euNQFa9erVU+0fP368TJgwQaZMmWKeU6RIEWnRooUcP348y9oKAADg6uDtxIkT8tBDD8nMmTPlqquu8u3XrJtm4IYMGSLt2rWTqlWryrx58+TPP/+UBQsWZGmbAQAAXBu89erVS1q3bi3NmzdPtX/Xrl1y4MABadmypW9f7ty5pVGjRrJmzZosaCkAAIDLx7y98847smHDBtMlmpYGbqpw4cKp9uv93bt3X/KYp0+fNpv/GAAAAIBIkWWZt8TERDM54Z///KdER0df8nnZsmVLdV+7U9PuSzvxIS4uzreVLFkyU9sNAADgyuBt/fr1cvDgQbnxxhslR44cZlu9erW88sor5rY34+bNwHnpz6TNxvkbPHiwHDt2zLdpkAgAABApsqzbtFmzZrJ5c+qp0I8++qhUrFhRnn76abn++uvN7NJly5ZJrVq1zONnzpwxAd6LL754yePquDjdAAAAIlGWBW/58+c3M0j95c2bV66++mrf/n79+smYMWOkfPnyZtPbefLkkY4dO2ZRqwEAALKWY4r0XsygQYMkOTlZevbs6SvSu3TpUhP4BSp+QbxExUQF9DNuKOzqhsKVWcWWgqlcj8uj2K6z3ndcD8BhwduqVatS3deJCSNGjDAbAAAAHFDnDQAAAOlH8AYAAGARgjcAAACLELwBAABYhOANAADAIgRvAAAAFiF4AwAAsIij6ryFUkLHBImNjRWno3Cls9hUwNam906w3HA93HJNKLYLBI/MGwAAgEUI3gAAACxC8AYAAGARgjcAAACLELwBAABYhOANAADAIgRvAAAAFiF4AwAAsIhrivSGG4UrnYXr4SxcD2exqaCwTW0N93mmJKeEpC1wHjJvAAAAFiF4AwAAsAjBGwAAgEUI3gAAACxC8AYAAGARgjcAAACLELwBAABYhOANAADAIhTpjZBikDa1NSMo7uosXA9n4Xo4S7g/l5OSkiSuR1zQrwl7kHkDAACwCMEbAACARQjeAAAALELwBgAAYBGCNwAAAIsQvAEAAFiE4A0AAMAirqnzFr8gXqJiohxf/8ymOk1uqElnU1vdcI42tdUtnwM2tTVYbjhH2IXMGwAAgEUI3gAAACxC8AYAAGARgjcAAACLELwBAABYhOANAADAIgRvAAAAFiF4AwAAsIhrivQmdEyQ2NjYsL2eG4o62nSONrXVDedoU1vdUFDYpra65X0XTFtTklPELbJly3bZxzt37ixz586VSOWa4A0AAESG/fv3+26/++67MmzYMNmxY4dvX0xMTKrnnz17VnLmzCmRgm5TAABglSJFivi2uLg4k4nz3j916pQUKFBA3nvvPWncuLFER0fLP//5TxkxYoTUrFkz1XEmTZok1113Xap9c+bMkUqVKpmfq1ixokybNk2chuANAABEnKefflr69Okj27dvl9tuuy1dPzNz5kwZMmSIjB492vzcmDFjZOjQoTJv3jxxErpNAQCAYyQlJaW6nzt3brMFql+/ftKuXbuAfmbUqFHy8ssv+36uTJkysm3bNpk+fboZR+cUZN4AAIBjlCxZ0nSFerexY8cGdZw6deoE9PxDhw5JYmKidO3aVfLly+fbXnjhBfn555/FSci8AQAAx9AAyr86RDBZN5U3b17xlz17dvF4PBdMZPA6f/68r+u0Xr16qZ4XFRUlTkLwBgAAHEMDt1CU9ipUqJAcOHDABHDeUiPffvut7/HChQtL8eLFZefOnfLQQw+JkxG8AQCAiNe4cWPTNTp+/Hhp3769LFmyRBYvXpwqUNQZqTrJQfe1atVKTp8+LevWrZMjR47IgAEDxCkI3i6DwpWRUbhSUUg0NLgeznrPcT3cXfxYB/rH9YgL+jUjXaVKlUzZD51BqhMT7r33Xhk4cKDMmDHD95xu3bpJnjx55KWXXpJBgwaZrtdq1aqZyQ9OQvAGAACs9cgjj5jNS+u2pR3b5tW9e3ez+Xv22WdT3e/YsaPZnIzZpgAAABYheAMAALAIwRsAAIBFCN4AAAAsQvAGAABgEYI3AAAAixC8AQAAWMQ1dd7iF8RLVExga5O5obCrGwpXZhVbCqZyPS6PYrvOet9xPQAybwAAAJGXeStYsGBAB9UFXzds2CClS5cOtl0AAAAINng7evSoTJo0SeLirrxmmi5J0bNnT0lJSUnPoQEAABCKMW8PPvigXHvttel67hNPPBFIGwAAAJCZwdv58+clEMePHw/o+QAAAEgfSoUAAABEcvA2b948+fTTT333Bw0aJAUKFJAGDRrI7t27M7t9AAAAyEjwNmbMGImJiTG3//vf/8qUKVNk/Pjxcs0110j//v0DPRwAAABCWaQ3MTFRypUrZ25/9NFH0r59e/n73/8uDRs2lMaNG4tTJXRMkNjYWHE6Clc6i00FbG167wTLDdfDLdeEYrtAGDNv+fLlk99//93cXrp0qTRv3tzcjo6OluTk5Aw0BQAAAJmeeWvRooV069ZNatWqJT/88IO0bt3a7N+6datcd911gR4OAAAAocy8TZ06VerXry+HDh2SRYsWydVXX232r1+/Xjp06BDo4QAAABDKzJvOLNVJCmmNHDky0EMBAAAg1Jk37Rp9/vnnzcQFAAAAODx4e/LJJ+Xjjz+WMmXKmPFv77zzjpw+fTo0rQMAAEDGgjddt1THt+lWuXJl6dOnjxQtWlR69+4tGzZsCPRwAAAACMfyWDVq1JDJkyfLr7/+KsOHD5c33nhD6tata/bPnj1bPB5PsIcGAABAZk1Y8Dp79qx8+OGHMmfOHFm2bJnEx8dL165dZd++fTJkyBBZvny5LFiw4LLHeO2118z2yy+/mPtVqlSRYcOGSatWrcx9DQB1IsSMGTPkyJEjUq9ePTPbVZ/ndBSudBauh7NwPZzFpoLCNrU13OeZkpwSkrYgAoI37RrVgO3tt9+WqKgo6dSpk0ycOFEqVqzoe07Lli3l1ltvveKxSpQoIePGjfOt2KDrpt59992yceNGE6DpslsTJkyQuXPnSoUKFeSFF14w4+x27Ngh+fPnD7TpAAAA7gvetGtUAyjNmN1zzz2SM2fOC56jY+EefPDBKx7rrrvuSnV/9OjR5rgJCQnmGJMmTTJZvHbt2vmCu8KFC5uM3uOPPx5o0wEAANwXvO3cuVNKly592efkzZvXZOcCkZKSIgsXLpSTJ0+aIsC7du2SAwcOmCyeV+7cuaVRo0ayZs2aSwZvOvPVf/ZrUlJSQO0AAACIqAkLVwrcArV582azXqoGZt27dzfj6DTrpoGb0kybP73vfexixo4dK3Fxcb6tZMmSmdpeAAAAq4I3zZD94x//kJtuukmKFCkiBQsWTLUF6oYbbpBvv/3WdJX26NFDOnfuLNu2bfM9ni1btlTP10kMaff5Gzx4sBw7dsy3UUwYAAC4OnjT2Z86ieD+++83wdGAAQPMmLTs2bPLiBEjAm5Arly5zISFOnXqmKyZtwSJBoYqbZbt4MGDF2Tj/GkGLzY2NtUGAADg2uBt/vz5MnPmTBk4cKDkyJHDLEavNd60xIdmzzJKM2s6Zk1XcNAATsuQeJ05c0ZWr14tDRo0yPDrAAAAuGLCgmbCqlX7q/6MjlXT7Ju68847ZejQoQEd69lnnzU13XRc2vHjx81SW6tWrZIlS5aYrtF+/frJmDFjpHz58mbT23ny5JGOHTtKONhUT8imtmYE9cGchevhLFwPZwn357JO0IvrERf0ayKCgzetzbZ//34pVaqU6e5cunSp1K5dW9auXWu6LAPx22+/mTpxejydXFC9enUTuGkpEjVo0CBJTk6Wnj17+or06utR4w0AALhVwMFb27ZtZcWKFSaQ6tu3r+k2nTVrluzZs0f69+8f0LH05y5Hs286ji6YsXQAAACRKODgTVdE8Grfvr3JxGndNc3CtWnTJrPbBwAAgMxY29RL1zTVDQAAAA4J3j755JN0H5DsGwAAQBYHb7qGadqxaFrSI+0+bxFfAAAAZGGdt/Pnz/s2ne1Zs2ZNWbx4sRw9etSUCtHbOuNUZ4oCAADAQWPetPba66+/LjfffLNv32233Wbqr/3973+X7du3Z3YbAQAAEGzw9vPPP5uabGnpvl9++UWcKn5BvETFRDm+eK1NRTbdUFDYpra64RxtaqtbPgdsamuw3HCOiPDlserWrWuyb1pY13/VhSeffNIsVg8AAAAHBW+zZ882i8OXLl3a1HbTTVdb0GDuSkV3AQAA3CA5OVn27dvnjG5TDdY2bdpkFoz//vvvzazTypUrS/PmzX0zTgEAANxq+fLlZkWqP//8U1q2bCkLFy4068FnWeZNaZCmjenTp49ZIkvXIiVwAwAAEBk4cKA8/PDDsm7dOlOpo1WrVpKUlBTe4O2VV16RU6dOpfugOhv1+PHjGWkXAACAlX788UcZNGiQ1KpVS/71r3+ZihzFixeX66+/XrZu3Wp6K/V2SIM3XXA+kGBMG3zo0KGgGwUAAGCr4sWL+8a75cqVSz777DOZP3++mfB57bXXmi7Vzp07h3bMm45ra9asmeTIkSPdg/QAAADc6L777pPp06dLw4YNzf2oqKhUy4f26tUrQ8dPVzQ2fPjwgA569913S8GCBYNtEwAAgLWGDRtmKnOESkiCNydK6JggsbGxYXs9NxR1tOkcbWqrG87Rpra6oaCwTW11y/sumLamJLO2uFPkzp1bSpYsGbLjB1wqBAAAAOmzd+9e+eSTT2TPnj1y5syZVI9NmDBBgkHwBgAAEAIrVqwwY93KlCkjO3bskKpVq5qlRHUuQe3atcNb5w0AAACXN3jwYLN86JYtWyQ6OloWLVokiYmJ0qhRIzOpIVgEbwAAACGwfft2X0kQrdih1Th0pYXnn39eXnzxxfAHb9pvqynAc+fOBf3iAAAAkSpv3rxy+vRpc7tYsWLy888/+x47fPhw+II3Xaera9euplpwlSpVzAA8pUtljRs3LuiGAAAARJL4+Hj5+uuvze3WrVubLtTRo0dLly5dzGNhC960//a7776TVatWmf5bL13q4d133w26IQAAAJFkwoQJUq9ePXN7xIgRZi14jZVKly4ts2bNCvq4Ac82/eijj8wLa8Tovxh95cqVU6UDAQAA3Ox6v/VLtcdy2rRpmXLcgIM3XbNU1+VK6+TJk6mCuUhA4crIKFypKCQaGlwPZ73nuB7uLn6clJQkcT3ign5NhCZ4W7t2rVx99dWp9h89etSUCtm5c2d4uk3r1q0rn376qe++N2CbOXOm1K9fP6hGAAAARJpffvlFUlIuXPlCJzH8+uuv4cu8jR07Vm6//XbZtm2bmWk6efJk2bp1q/z3v/+V1atXB90QAACASPDJJ5/4bn/++ecSF/f/GVEN5rR473XXXRe+4K1BgwZm5sQ//vEPKVu2rCxdutSk/jR4q1Yt+BQxAABAJLjnnnt8vZPeOm9eOXPmNIHbyy+/HPTxg1oeS4O0efPmBf2iAAAAker8+fPm/7oslo55u+aaazL1+EGvbXrw4EGzeRvoVb169cxoFwAAgNV27doVkuMGHLytX7/epAB1yQddWNWfpgcvNjAPAADAjU6ePGnmBOiiBro6lT9d4CAswdujjz4qFSpUMMXlChcuHHHlQQAAADLDxo0b5Y477jCrU2kQV7BgQbMsltZ807JrYQveNAX4wQcfSLly5YJ6QQAAADfo37+/3HXXXfLaa69JgQIFJCEhwUxYePjhh6Vv375BHzfg4K1Zs2ZmeSzbgrf4BfESFRMV0M+4obCrGwpXZhVbCqZyPS6PYrvOet9xPWCTb7/9VqZPny5RUVFm0/puWrh3/PjxZghau3btwhO8vfHGG+YFt2zZIlWrVjURpL82bdoE1RAAAIBIkjNnTt/wMh1qpuPeKlWqZOq+6e1gBRy8rVmzRr766itZvHjxBY8xYQEAAOAvtWrVknXr1pm5Ak2aNJFhw4aZMW9vvfVWhmrjBrw8lg6u69Spk+zfv9+UCfHfmGkKAADwlzFjxkjRokXN7VGjRpk1Tnv06GFKrc2YMUPClnn7/fffzQA8Tf8BAADg4urUqeO7XahQIfnss88kMwScedPBdStXrsyUFwcAAECIM2/abzt48GAz7k37a9NOWAi2ZgkAAEAkjHPLls4auBs2bAjfbNN8+fKZasG6+dPGErwBAAC3L0qvTp06JdOmTZPKlStL/fr1zT6t9bZ161bp2bOnBCuoIr0AAAC40PDhw323u3XrZpJaOlkh7XMSExMl7AvT2yahY4LExsaK01G40llsKmBr03snWG64Hm65JhTbhRssXLjQlApJS1dY0MkMs2fPDl3wNmDAABM15s2b19y+nAkTJgTVEAAAgEgSExNj5giUL18+1X7dFx0dHfRxc6R3YdWzZ8/6bgMAAODy+vXrZ+q6rV+/XuLj431j3jTjpgV7Qxq8+ZcGoUwIAADAlT3zzDNmLdPJkyfLggULzD5dHmvu3Lly//33S9jqvHXp0kWOHz9+wf6TJ0+axwAAAPAXDdK+/vpr+eOPP8ymtzMSuAUVvM2bN0+Sk5Mv2K/73nzzzQw1BgAAAJk02zQpKUk8Ho/ZNPPmP9BO1zTVJR+uvfba9B4OAAAAoQzeChQoYIrw6qarLKSl+0eOHBlMGwAAAJDZwZtOVNCsW9OmTWXRokVSsGBB32O5cuWS0qVLS7FixdJ7OAAAAIQyeGvUqJFvhYVSpUqle90ut6JwpbNwPZyF6+EsNhUUtqmt4T7PlOSUkLQFmUeHmW3evNkkvK666qrwTVjQFyRwAwAAuHKdt1mzZvkCN02E1a5dW0qWLCmrVq2SsAVvAAAAuLL3339fatSoYW7/61//Mr2X33//vQnqhgwZIsEieAMAAAiBw4cPS5EiRcxtrcpx3333mUmfXbt2Nd2nwSJ4AwAACIHChQvLtm3bTJfpkiVLpHnz5mb/n3/+KVFRUaGfsAAAAID0e/TRR81qCkWLFjXzBVq0aGH2f/PNN1KxYkUJafBWq1atdE9S2LBhQ9CNAQAAiBQjRoyQqlWrSmJioukyzZ07t9mvWTdd9zSkwds999zju33q1CmZNm2aVK5cWerXr2/2JSQkyNatW6Vnz55BNwQAACDStG/f/oJ9nTt3ztAxs3m08m4AunXrZtJ/o0aNSrV/+PDhJrKcPXu2OIku6xUXFyfHjh2T2NjYiK0nZFNbM4L6YM7C9XAWroezhPtzOSPfd8HI7NfzHk8H9GtmqlevXmaz3YoVK8x28OBBOX/+fKrHgo2ZAh7ztnDhQlm3bt0F+x9++GGpU6eO44I3AABgj7Vr14Yl+AwHXTb0+eefN/GRd9xbZgg4eIuJiZGvvvpKypcvn2q/7vNfrB4AAMDNXn/9dZk7d6506tQpU48bcPCmheV69Ogh69evl/j4eN+YN824DRs2LFMbBwAAYKszZ85IgwYNMv24AQdvOjvi+uuvl8mTJ8uCBQvMvkqVKpnIUqfDAgAAQMw8AY2Vhg4dmqm/jqDqvGmQRqAGAABwaVqhY8aMGbJ8+XKpXr265MyZM9XjEyZMkGBQpBcAACAENm3aJDVr1jS3t2zZkuqxjExeCDh40yUeJk6cKO+9957s2bPH9Of6++OPP4JuDAAAQKRYuXJlSI6bPZhpr5rm025Tre0yYMAAadeunWTPnt1UEgYAAEDoBJx5mz9/vsycOVNat25tArkOHTpI2bJlTV+uzjrt06ePOFH8gniJiolyfPFam4psuqGgsE1tdcM52tRWt3wO2NTWYLnhHJF5NKGlkzi1Vp3evpwPPvggPMHbgQMHpFq1v97I+fLlM9k3deedd2b6bAoAAACbxMXF+caz6e1QCDh4K1GihOzfv19KlSol5cqVk6VLl0rt2rVNRWTvgqsAAABuNGfOnIveztIxb23btjVrdKm+ffuabJuutvC3v/1NunTpEoo2AgAAINjM27hx43y327dvbzJxa9asMVm4Nm3aBHo4AACAiPX+++9fskLHhg0bwpN5S0uXyNIZpwRuAAAA/++VV16RRx99VK699lrZuHGj3HTTTXL11VfLzp07pVWrVhKsoIK3t956Sxo2bCjFihWT3bt3m32TJk2Sjz/+OOiGAAAARJJp06aZFRamTJkiuXLlkkGDBsmyZctMZQ7vhM+wBG+vvfaaybTdcccdcvToUVO0VxUoUMAEcAAAABDTVepdmD4mJkaOHz9ubnfq1Enefvvt8AVvr776qqnzNmTIEImK+v+6aXXq1JHNm6lpAwAAoIoUKSK///67uV26dGlTD1ft2rVLPB6PhG3Cgr5grVq1LtivZUJOnjwZ0LHGjh1rCtR9//33JiLV6PTFF1+UG264wfccPTktBqxpxyNHjki9evVk6tSpUqVKlYBeK6FjgimYFy5uKOpo0zna1FY3nKNNbXVDQWGb2uqW910wbU1J/qsnDM7RtGlT+de//mVKqnXt2lX69+9vJjCsW7fuigV8MzXzVqZMGfn2228v2L948WKpXLlyQMdavXq19OrVy0Si2gd87tw5admyZaogcPz48WY5Lu0v1lpyGsW2aNHCl3oEAABwIk08aU+l6t69u1l5oVKlSiYppcPQwpZ5e+qpp0zAderUKZMV+9///mf6bTWL9sYbbwR0rCVLlqS6r8XsdEbG+vXr5dZbbzXH13F0euLeCHXevHlSuHBhWbBggTz++OOBNh8AACAs9u7dKyVLlvTd13XhddP4JjEx0Sx4EJbgTae8aoZMZ0z8+eef0rFjRylevLhMnjxZHnzwQckI78yLggUL+rpodTkuzcb5d882atTI1Ja7WPB2+vRps3klJSVlqE0AAADB0N5KXZVKE1P+/vjjD/OYd9JnWEqFPPbYY6ZEyMGDB01wpdGj9uVmhEahOov15ptvlqpVq5p9emylmTZ/et/7WFqaAdS1xLybf8QLAAAQLhrbeNc59XfixAmJjo4O+rgBZ978XXPNNZJZevfuLZs2bZKvvvrqgsfSnvilfhlq8ODBJgj0z7wRwAEAgHDxxiEaq+gyonny5PE9ptm2b775RmrWrBna4E1nl14qWEormKUennjiCfnkk0/kyy+/NMtteenkBKVZtqJFi/r2a8YvbTbOv1tVNwAAgKygqyl4k01aRk0L9Hrp7Ro1asjAgQNDG7zdc889Egp6Uhq4ffjhh7Jq1SrT/+tP72sApzNRveVJdF0wnaWqJUUAAACcZuXKlb55AjonILNLlaUreBs+fLiEgs5a1VmjuqxW/vz5fePYdKya1n3TbF+/fv1kzJgxUr58ebPpbU0/6kQJAAAAp9IqGqGQoTFvGeWtcdK4ceMLTvaRRx4xt3VWa3JysvTs2dNXpHfp0qUm2As1CldGRuFKRSHR0OB6OOs9x/Vwd/FjHeMd1yMu6NdE5tO6tePGjZMVK1aYIV/nz59P9bguUG9d8JaepSE0+zZixAizAQAA2KJbt25mqJeuZapj99M7f8DRwRsAAECkWrx4sXz66afSsGHDTD1uUHXeAAAAcHlXXXWVb+GBLA3enn/+ebOyQlo6Lk0fAwAAgMioUaNk2LBhF42bwhq86WKqWhk4LW2YPgYAAACRl19+WT7//HNTm7ZatWpSu3btVFvYxrxdanWD7777LiSpQQAAABvdE6I6uTkC6bfVoE23ChUqpArgdKkHzcZ17949JI0EAACwzfAQ1clNd/A2adIkk3Xr0qWL6R7VQrr+Sz1cd911Ur9+/ZA0EgAAwFbr16+X7du3m8RX5cqVfatGhTx469y5s2/JKp3ymiOHXVVG4hfES1RMVEA/44bCrm4oXJlVbCmYyvW4PIrtOut9x/WATbQw74MPPmiWAC1QoIBJgh07dkyaNGki77zzjhQqVCg8ExZ0ZQONHr10aSvt03322WfNuqMAAAAQs367rnyxdetW+eOPP8xKUVu2bDH7+vTpE/SvKODg7fHHH5cffvjBt6zDAw88YNYaXbhwoVnKCgAAACJLliwxS4FWqlTJ9+vQbtOpU6eaAr5hC940cKtZs6a5rQFbo0aNzOLyc+fOlUWLFnGtAAAARMxapjlz5rzgd6H70q5zGtLgTftrvS+4fPlyueOOO8ztkiVLyuHDh7lYAAAAItK0aVPp27ev7Nu3z/f7+PXXX6V///7SrFmz8AVvderUkRdeeEHeeusts9hq69atzf5du3aZInQAAAAQmTJlihw/ftxU5ChbtqyUK1fOTPzUfa+++mrQv6KAp4xqyZCHHnpIPvroIxkyZIhpiHr//felQYMGXCsAAAD5q1dyw4YNsmzZMvn+++9N76WOeWvevHmGfj8BB2/Vq1eXzZsvLGXw0ksvSVRUYKU4AAAAIl2LFi3MllkC7ja9lOjo6IsOygMAAHCTL774wmTYtCRIWlrnrUqVKvKf//wnfJk3XQpr4sSJ8t5778mePXsuqO2mdUycKKFjgsTGxorTUbjSWWwqYGvTeydYbrgebrkmFNtFJJs0aZI89thjF407dIUqLbs2YcIEueWWW8KTedOlsfQF77//fhM9DhgwQNq1ayfZs2eXESNGBNUIAACASPHdd9/J7bfffsnHW7ZsaZbMClbAwdv8+fNl5syZMnDgQLNEVocOHeSNN96QYcOGSUJCQtANAQAAiAS//fbbZYeSafx06NCh8AVvBw4ckGrV/kqT58uXz2Tf1J133imffvpp0A0BAACIBMWLF7/o5E6vTZs2SdGiRcMXvJUoUUL2799vbmuZkKVLl5rba9euldy5cwfdEAAAgEhwxx13mB7JU6dOXfBYcnKyDB8+3CS9wjZhoW3btrJixQqpV6+eqRqs3aazZs0ykxe0YjAAAICbPffcc/LBBx9IhQoVpHfv3nLDDTdItmzZZPv27WZdU538qbVywxa8jRs3zne7ffv2pgDd119/bbJwbdq0CbohAAAAkaBw4cKyZs0a6dGjhwwePNgU51UawN12220ybdq0DK1KFXDwlpZm4HQDAADAX0qXLi2fffaZHDlyRH766ScTwJUvX16uuuoqyaiAg7exY8eaaLFLly6p9s+ePdvMnHj66acz3KhIQO0jZ+F6OAvXw1lsqklnU1vDfZ4pySkhaQsyRoO1unXrSmYKeMLC9OnTpWLFihfs12rBr7/+ema1CwAAAJlVKuRi01sLFSrkm4UKAAAAhwRv3gkKaem+YsWKZVa7AAAAkBlj3rp16yb9+vWTs2fPStOmTc0+LR0yaNAgefLJJwM9HAAAAEIZvGmQpovP9+zZ07cofXR0tJmooNNhAQAA4KDgTWuUvPjiizJ06FBTbC4mJsZMfWV1BQAAAAeOefOfuKAZuLJly5rAzVuADgAAAA4K3n7//Xdp1qyZWfJB1+7yzjDVsXCMeQMAAHBYt6muX5ozZ06zlmmlSpV8+x944AHz2MsvvyyRwqZikDa1NSMo7uosXA9n4Xo4S7g/l5OSkiSuR1zQr4kIDt6WLl0qn3/+uZQoUSLVfh33tnv37sxsGwAAADLabXry5EnJkyfPBfsPHz7MpAUAAACnBW+33nqrvPnmm6lmn54/f15eeukladKkSWa3DwAAABnpNtUgrXHjxrJu3TpT503rvm3dutXMPL3YygsAAADIwsxb5cqVZdOmTXLTTTdJixYtTDdqu3btZOPGjaZsCAAAABySedMlsVq2bCnTp0+XkSNHhq5VAAAAyHjmTUuEbNmyxYxzAwAAgAXdpn/7299k1qxZoWkNAAAAMnfCgk5SeOONN2TZsmVSp04dyZs3b6rHJ0yYIE4UvyBeomKiHF+81qYim24oKGxTW91wjja11S2fAza1NVhuOEdEePCm3aa1a9c2t3/44YdUj9GdCgAA4LDgbeXKlaFpCQAAADJ/zJu/vXv3yq+//pqRQwAAACCUwZuupvD8889LXFyclC5dWkqVKiUFChSQUaNGmccAAADgoG7TIUOGmNmm48aNk4YNG4rH4zErK4wYMUJOnTolo0ePDk1LAQAAEHjwNm/ePDPbtE2bNr59NWrUkOLFi0vPnj0J3gAAAJzUbaprmFasWPGC/bpPHwMAAICDgjfNsk2ZMuWC/bpPHwMAAICDuk3Hjx8vrVu3luXLl0v9+vVNbbc1a9ZIYmKifPbZZ+JUCR0TJDY2Nmyv54aijjado01tdcM52tRWNxQUtqmtbnnfBdPWlOSUkLQFEZB5a9SokSnO27ZtWzl69KjpKm3Xrp3s2LFDbrnlltC0EgAAAIFl3nbu3CllypQxmbZixYoxMQEAAMDJmbfy5cvLoUOHfPcfeOAB+e2330LVLgAAAGQkeNN6bv50fNvJkyfT++MAAADI6uWxAAAA4NDgTce66ZZ2HwAAABw4YUG7TR955BHJnTu3ua9LYXXv3l3y5s2b6nkffPBB5rcSAAAAgQVvnTt3TnX/4YcfTu+PAgAAINzB25w5c8RtKFwZGYUrFYVEQ4Pr4az3HNfD3cWPk5KSJK5HXNCvCXswYQEAAMAiBG8AAAAWIXgDAACwCMEbAACARQjeAAAALELwBgAAYBGCNwAAAIsQvAEAANeYO3euFChQQFxRpNd28QviJSomKqCfcUNhVzcUrswqthRM5XpcHsV2nfW+43rAS5fsnDdvnqT1448/Srly5SSSuSZ4AwAAkeX222+/YAWoQoUKSaSj2xQAAFgpd+7cUqRIkVTb5MmTpVq1apI3b14pWbKk9OzZU06cOHHJY3z33XfSpEkTyZ8/v8TGxsqNN94o69at8z2+Zs0aufXWWyUmJsYcr0+fPnLy5EnJSgRvAADAMXSNVv/t9OnTAf189uzZ5ZVXXpEtW7aYbtUvvvhCBg0adMnnP/TQQ1KiRAlZu3atrF+/Xp555hnJmTOneWzz5s1y2223Sbt27WTTpk3y7rvvyldffSW9e/eWrES3KQAAcAzNbvkbPny4jBgx4qLP/fe//y358uXz3W/VqpUsXLjQd79MmTIyatQo6dGjh0ybNu2ix9izZ4889dRTUrFiRXO/fPnyvsdeeukl6dixo/Tr18/3mAaGjRo1ktdee02io6MlKxC8AQAAx0hMTDTdl/5do5fSpEkTE0R5aVfpypUrZcyYMbJt2zaTuTt37pycOnXKdHXq42kNGDBAunXrJm+99ZY0b95c7rvvPilbtqx5TDNxP/30k8yfP9/3fI/HI+fPn5ddu3ZJpUqVJCvQbQoAABxDAzf/7XLBW968ec3MUu925swZueOOO6Rq1aqyaNEiE3xNnTrVPPfs2bMXPYZm9bZu3SqtW7c2XayVK1eWDz/80DymQdrjjz8u3377rW/TMXI6o9Ub4GUFMm8AACAirFu3zmTaXn75ZTP2Tb333ntX/LkKFSqYrX///tKhQwczg7Vt27ZSu3ZtE9g5rfRIlgZvX375pelP1sh4//79JtK95557UqUmR44cKTNmzJAjR45IvXr1TARdpUqVgF8roWNCqjSsU1H7yFlsqoFm03snWG64Hm65JtRrQyiULVvWBG+vvvqq3HXXXfL111/L66+/fsnnJycnm/Fu7du3N+Pj9u7dayYu3Hvvvebxp59+WuLj46VXr17y2GOPmUzf9u3bZdmyZeY1skqWdptq/3ONGjVkypQpF318/PjxMmHCBPO4/jJ1CnCLFi3k+PHjYW8rAABwtpo1a5q44cUXXzRdpzpWbezYsZd8flRUlPz+++/yt7/9zWTe7r//fjPpQRNHqnr16rJ69WrTTXrLLbdIrVq1ZOjQoVK0aFFxbeZNf0G6XYxm3SZNmiRDhgwxU3SVTvktXLiwLFiwwPRBAwAA9y5zdTHa9ambv06dOqVamUE3lStXLnn77bflcurWrStLly4VJ3HshAWdxXHgwAFp2bKlb58OWtTpuVowDwAAwI0cO2FBAzelmTZ/en/37t2X/Dkt5udf0E+nCQMAAEQKx2bevLJly3ZBd2raff60bzsuLs63pS32BwAAYDPHBm86OcE/A+d18ODBC7Jx/gYPHizHjh3zbVrsDwAAIFI4NnjTKbsawOl0XC8tvqezPho0aHDJn9NxcWkL/AEAAESKLB3zduLECbPshP8kBa1eXLBgQSlVqpRZS0yXuNC1xHTT23ny5DHrjAEAALhRjqyuhKzrkvmvL6Y6d+5spgAPGjTIFNDr2bOnr0ivTtfNnz+/OB2FK52F6+EsXA9nsamgsE1tDfd5piSnhKQtcJ4sDd4aN25sJiBcik5M0DXHdAMAAICDx7wBAADgQgRvAAAAFiF4AwAAsAjBGwAAgEUI3gAAACxC8AYAAGARgjcAAACLZGmdN6ezqRikTW3NCIq7OgvXw1m4Hs4S7s/lpKQkiesRF/Rrwh5k3gAAACxC8AYAAGARgjcAAACLELwBAABYhOANAADAIgRvAAAAFiF4AwAAsAjBGwAAgEVcU6Q3fkG8RMVEOb54rU1FNt1QUNimtrrhHG1qq1s+B2xqa7DccI6wC5k3AAAAixC8AQAAWITgDQAAwCIEbwAAABYheAMAALAIwRsAAIBFCN4AAAAsQvAGAABgEdcU6U3omCCxsbFhez03FHW06RxtaqsbztGmtrqhoLBNbXXL+y6YtqYkp4SkLXAeMm8AAAAWIXgDAACwCMEbAACARQjeAAAALELwBgAAYBGCNwAAAIsQvAEAAFiE4A0AAMAirinSGwwKV0ZG4UpFIdHQ4Ho46z3H9XB38eOkpCSJ6xEX9GvCHmTeAAAALELwBgAAYBGCNwAAAIsQvAEAAFiE4A0AAMAiBG8AAAAWIXgDAACwCMEbAACARVxTpDd+QbxExUQF9DNuKOzqhsKVWcWWgqlcj8uj2K6z3ndcD4DMGwAAgFXoNgUAALAIwRsAAIBFCN4AAAAsQvAGAABgEYI3AAAAixC8AQAAWMQ1dd4SOiZIbGysOB21j5zFphpoNr13guWG6+GWa0K9NiB4ZN4AAAAsQvAGAABgEYI3AAAAixC8AQAAWITgDQAAwCIEbwAAABYheAMAALAIwRsAAIBFXFOkN9woXOksXA9n4Xo4i00FhW1qa7jPMyU5JSRtgfOQeQMAALAIwRsAAIBFCN4AAAAsQvAGAABgEYI3AAAAixC8AQAAWITgDQAAwCIEbwAAABahSG+EFIO0qa0ZQXFXZ+F6OAvXw1nC/bmclJQkcT3ign5N2IPMGwAAgEUI3gAAACxC8AYAAGARgjcAAACLELwBAABYhOANAADAIgRvAAAAFiF4AwAAsIgVRXqnTZsmL730kuzfv1+qVKkikyZNkltuuSWgY8QviJeomCjHF6+1qcimGwoK29RWN5yjTW11y+eATW0NlhvOEXZxfObt3XfflX79+smQIUNk48aNJmhr1aqV7NmzJ6ubBgAAEHaOD94mTJggXbt2lW7dukmlSpVM1q1kyZLy2muvZXXTAAAAws7RwduZM2dk/fr10rJly1T79f6aNWsu+jOnT58267v5bwAAAJHC0cHb4cOHJSUlRQoXLpxqv94/cODARX9m7NixEhcX59s0SwcAABApHB28eWXLli3VfY/Hc8E+r8GDB8uxY8d8W2JiYphaCQAA4PLZptdcc41ERUVdkGU7ePDgBdk4r9y5c5sNAAAgEjk685YrVy658cYbZdmyZan26/0GDRpkWbsAAACyiqMzb2rAgAHSqVMnqVOnjtSvX19mzJhhyoR07949XT+vXawqJTkl4NfOiskOwbRT2TQxI9hzzIrztKmtbjhHm9rqls8Bm9oa6efofT3v9164X9cpx3GDbJ5wX+Ugi/SOHz/eFOmtWrWqTJw4UW699dZ0/ezevXuZtAAAcA0d612iRImQv86pU6ekTJkyl5xAGIzY2FgpWrSoZM+eXXr16mU2WBq8ZcT58+dl3759kj9//gsmOWiUr7NR9Y2ub5hI5IZzdMt5co6Rg2sZGZx2HfXr/Pjx41KsWDET/ISDBnBa1iszh0tFR0dn2vEileO7TTNK38BX+gtE/9E54R9eKLnhHN1ynpxj5OBaRgYnXUctkRVOGmgRbIWfoycsAAAAIDWCNwAAAIu4OnjTenDDhw+P6LpwbjhHt5wn5xg5uJaRwQ3XEc4U8RMWAAAAIomrM28AAAC2IXgDAACwCMEbAACARSI6eHvkkUfknnvuEbfQlSi02rXW3NE1Yf/zn/+k+l1okWL/LT4+XiLtPH/77TdzrlqkMk+ePHL77bfLjz/+KDb58ssv5a677jLnoNfpo48+SvX4iBEjpGLFipI3b1656qqrpHnz5vLNN99IJJ1j2veqd3vppZfEFmPHjpW6deuaAuHXXnut+SzasWNHqud88MEHctttt8k111xjzu/bb78Vm6TnHG3/7EnPOUbC5w7sEtHBm5u8++670q9fPxkyZIhs3LhRbrnlFmnVqpVZB9ZLP1B0iTHv9tlnn0kknafOvdEP1p07d8rHH39sHi9durQJbk6ePCm20LbWqFFDpkyZctHHK1SoYB7bvHmzfPXVV3LddddJy5Yt5dChQxIp5+j/PtVt9uzZ5kv/3nvvFVusXr3aLO2TkJAgy5Ytk3Pnzpnr5P9e1NsNGzaUcePGiY3Sc462f/Zc6Rwj5XMHlvFEsM6dO3vuvvtuc3vx4sWehg0beuLi4jwFCxb0tG7d2vPTTz/5nrtr1y6ddetZtGiRp3Hjxp6YmBhP9erVPWvWrPHY4KabbvJ079491b6KFSt6nnnmmQt+Fza73Hnu2LHDXMMtW7b4Hjt37py53jNnzvTYSM/nww8/vOxzjh07Zp63fPlyT6Seo753mzZt6rHZwYMHzbmuXr36gse8nz8bN270RNo5Rspnz6XOMRI/d+B8rsm86V9AAwYMkLVr18qKFSvMsllt27Y1a5/604zOwIEDTfeFZjg6dOhg/tJyMl1Xbv369eavQX96f82aNb77q1atMml/Pa/HHntMDh48KDa50nmePn3a3PdfqiUqKsqslacZqkikv5MZM2aYJXE0kxWJtEvq008/la5du4rNjh07Zv5fsGBBiVSXOkfbP3sud45u/NxB1nNN8KbdLe3atZPy5ctLzZo1ZdasWabbadu2bamep4Fb69atzYfMyJEjZffu3fLTTz+Jkx0+fFhSUlKkcOHCqfbr/QMHDpjb2rU4f/58+eKLL+Tll182QWzTpk19Hzw2uNJ56jgw7a4YPHiwHDlyxAQ22h2lj2lXTST597//Lfny5TNfGBMnTjTdOTpuKhLNmzfPjDfSf7+20gSj/vF48803S9WqVSUSXeocI+Gz53Ln6KbPHThHxC9M7/Xzzz/L0KFDzbgFDQK8GTcdK+X/QVO9enXf7aJFi5r/61+J+g/U6XRMUNoPGu++Bx54wLdfz7dOnTrmA0czGrZ9KV7qPHPmzCmLFi0yGRr9q1j/+tVxJ/rlEWmaNGlissP6Xp45c6bcf//9ZtKCZjcijY53e+ihh6xe/Lp3796yadOmiM7EXOocI+mz52Ln6KbPHTiHazJvOrPt999/N190+iXnnZ2nfyX503+IaYOEtF2rTqMZF/3A8GbZvDToTJul8g9M9QPUphlR6TlPnX2qQc3Ro0fNX71Lliwx111np0YSnWlarlw5M2tPs8g5cuQw/480OpNYZ/Z169ZNbPXEE0/IJ598IitXrpQSJUpIJArkHG387LnSObrlcwfO4YrgTf8Rbd++XZ577jlp1qyZVKpUyaS3I4WOrdAPD+0686f3GzRocMnfSWJioi+7GGnnqWPAChUqZL4g1q1bJ3fffbdEMs0+2tgNdSUakOo1t3E8n14TzdRoORDtMozEL/JgztG2z55AztFtnzvIOq7oNtVaWFdffbUZ2K0fGNpV+swzz0gk0XEYnTp1Ml0S9evXN+eq59m9e3c5ceKEqQ2m4/70/H/55Rd59tlnTSZLJ21EynmqhQsXmg/PUqVKmTGNffv2NdP4005ycDK9Xv7jLHft2mX+qtcuGX0fjx49Wtq0aWOupX4Rat27vXv3yn333SeRcI567VRSUpK5njpOykZaXmLBggWmfISO2fNmjPULPiYmxtz+448/zPt337595r63fliRIkXMZvs5RsJnT3quYyR87sAyngjWqVMnz7333mtuL1u2zFOpUiVP7ty5TQmQVatWpSpRcLGp+keOHDH7Vq5c6bHB1KlTPaVLl/bkypXLU7t2bd9U9j///NPTsmVLT6FChTw5c+b0lCpVykzf37Nnj8dGlzpPNXnyZE+JEiV85/ncc895Tp8+7bGJvt/0fZd202uWnJzsadu2radYsWLm/IsWLepp06aN53//+58nUs7Ra/r06aZkz9GjRz02utj56TZnzhzfc/T2xZ4zfPhwTyScYyR89qTnOkbC5w7skk3/IxFKC0PquKBLFQIFAACwTUSOedPxbDqTSWsL6awfAACASBGRY966dOliagk9+eSTDBgFAAARJaK7TQEAACJNRHabAgAARCqCNwAAAIsQvAEAAFiE4A0AAMAiBG8AAAAWIXgD4DiPPPKIWV7Iq3HjxtKvX7/L/sx1110nkyZNCuh1dLmmbNmymaW5EFnGjh0rdevWNUtaXXvtteb95F1+zEuLLejyXcWKFTNLXen7bOvWrameo0vw6f7Y2FjzXtHF5y9WW1SX7dMls3TT2xd7Xlr62jVr1syEs4XbELwBIQo+9IM+7ea/nqdt5s6dKwUKFMiS19ZFwUeNGpWpAaEqWbKk7N+/X6pWrZrBFsJpVq9ebdYlTUhIkGXLlsm5c+fMWqMnT570PWf8+PEyYcIEswqP1gbV9WRbtGghx48f9z3nzz//NKv16Jqsl9KxY0fzB8CSJUvMprc1gANCheANCBH9wNfAwH8rU6ZMUMc6c+aMuJkuWK8ZlMwWFRVlvrBz5IjIeuWupkGUBuxVqlSRGjVqyJw5c2TPnj2yfv16X9ZNM7VDhgyRdu3amQB+3rx5JljThei9NOP7zDPPSHx8/EVfZ/v27ea13njjDalfv77ZZs6cKf/+978vyPRdiQaQGjxec801JoPXqFEj2bBhQ6rn6B+B+lpt27aVPHnySPny5eWTTz4J6ncEexG8ASGSO3duExj4bxoseLMCN910k3lO0aJFzZeDZga8tJumd+/eMmDAAPNBrh/oatu2bXLHHXdIvnz5pHDhwuav+8OHD/t+7vz58/Liiy+aNX312KVKlZLRo0f7Hn/66aelQoUK5kP/+uuvl6FDh8rZs2d9j3/33XfSpEkTEyhpN9GNN94o69atM0vNPfroo3Ls2DFfFlG7fNLSLyt97Pvvv0+1X7Mb2q2pX5gpKSnStWtXE8hqV9UNN9wgkydPvuzvMm236cGDB+Wuu+4yP6/HmT9//mV/XtuqX8wff/yxr/16Tmm7TXWf3v/888+lVq1a5vhNmzY1r7d48WKpVKmS+b106NDBfMl76XlpFkd/p/ozGiy8//77l20Twkvfu94/BNSuXbvkwIEDJhvnpf9mNGBas2ZNuo/73//+1wRa9erV8+3TQE/3BXIcpRm/zp07y3/+8x+TMdTATP+9+2cC1ciRI+X++++XTZs2mccfeugh+eOPPwJ6LdiNPzeBMPv111/NB65mBd58800T6Dz22GMSHR2dKiDSYKNHjx7y9ddfm+BAM3f6xaLP1WAoOTnZBGP6If7FF1+Ynxk8eLD5q3/ixIly8803m5/xD6Q0KNPuTx3js3nzZnMs3Tdo0CDzuH4JaNDy2muvmUBTg5qcOXNKgwYNTJZi2LBhvmyCBpBpaSCmAZ8GU/7dnJrJ0K4lDYw0wCxRooS89957JjDVL7i///3vJojVc0kP/d0lJiaa886VK5f06dPHBFiXMnDgQJMhSUpKMhkY75f4vn37Lvp8vQ7alaZBrrZJN/1i1/M4ceKEyXq8+uqr5vevnnvuOdO1q783/cL98ssv5eGHH5ZChQqZa4aspf9+9A8h/Tfh7SLXwE3pH0H+9P7u3bvTfWw9jo6pS0v3eV8jvfQPBX/Tp0+Xq666yvyxd+edd6Z6/+sfEGrMmDHmvfi///3PZPvhEro8FoDM1blzZ09UVJQnb968vq19+/bmsWeffdZzww03eM6fP+97/tSpUz358uXzpKSkmPuNGjXy1KxZM9Uxhw4d6mnZsmWqfYmJibq8nWfHjh2epKQkT+7cuT0zZ85MdzvHjx/vufHGG3338+fP75k7d+5FnztnzhxPXFzcFY85YcIEz/XXX++7r23TNm7duvWSP9OzZ0/Pvffem+r3d/fdd/vu6++jb9++qY6XkJDge3z79u1m38SJEy/5GmmPqXbt2mV+buPGjeb+ypUrzf3ly5f7njN27Fiz7+eff/bte/zxxz233XabuX3ixAlPdHS0Z82aNamO3bVrV0+HDh0u2R6Ej76/Spcubf69eH399dfmuu7bty/Vc7t16+a7tv68740jR46k2j969GhPhQoVLnh+uXLlzHtHVa5c2fc5cPvtt/ueM3z4cE+NGjV893/77Tfz3ipfvrwnNjbWPD9btmzm88FL2/Dee++lei197rx58wL8rcBmZN6AENHuR83EeOXNm9f8XzNAOi5Gs1BeDRs2NBmdvXv3mq5OVadOnVTH07E6K1euvGjG6+effzaz206fPi3NmjW7ZJu0K08zaDpxQl9Pu2q1G9BLsxPdunWTt956S5o3by733XeflC1bNqDzfvDBB+Wpp54y3T7afaRZOJ1RV7lyZd9zXn/9dTNuRzMcmkHUMX3pnXWnvz8do+b/+6lYsWKmTqaoXr16qkyMt5vZf59mOrxd2adOnfJ1bXvpOWkWE1nriSeeMGPCNBuqGV8vHcagNDumWV8vzeCmzcZdjh7nt99+u2D/oUOHfMf57LPPfMMTtFv9UjSjpj+n/0ZLly5tsr36WZF2zKtmw/15M9pwD8a8ASGiwZqOPfNu3i8I/ePZP3Dz7lP++73Bnpd+OOs4L+3K9N9+/PFHufXWWy/7paA0mNLAqlWrVmYw9caNG81gbf8vBu0u1FIJrVu3Nl2SGnB9+OGHAZ23nqcGrt5B32+//bbpQvTS7tL+/ftLly5dZOnSpeYcdDxdeidlXOx3ldn8vxz1dS73Zen9/6effprqumhQx7i3rKPvEx03qt3Z+l5OO1lI72vgpTNRvfQ9qF2UOkwgvTS40vF03mBeffPNN2af9zgaiHk/B4oXL37JY+lYNx0CoMMqdKKFBm/+Y1oBLzJvQJhpQLRo0aJUQZyO+9KxZ5f7YK9du7b5OR34f7HZkTrWSgO4FStWmOxZWjp2Tr9ENGDzutjYHp3QoJsGWDquRseI6RgvHVumkw3SQ8fO6Xgw/XnNCmrQ6P8FpV9qPXv29O3T56SXThrQjKFOpNBJH0rH4V2prlYg7Q/0euqXrM5kZHybc2iZEP0DQiep6L8t7/gznUig/070355OgtExY/pvRze9rVlWHZ/ppT+nm7fMj44V1eNphlzHTer7Ucea6fhRHaOmdAynjlHTMaCB0OBOs96aVdbxmZrBvtIfZXAnMm9AmGnQooPttTtHJxPol8vw4cNNl2X27Nkv+2WkM8o0INK/8nfu3GkyV5rB0qBEJzxowKSTD3QihAZEmm2bNWuW74tBA4x33nnHPPbKK6+kyqpp96VmKnTGpQZ1Guxp6QL9clIaNGpXqwaHmg3wn22ZlpZe0C8fnXChWTj/oFTboYGXzuj84YcfzIxXfZ300i9E75elZji0O1mD1St9yWn7dXaeBnrafv9ZthmhX+Q6IUKDXZ1kor9bzWpOnTrV3EfW0CELmv3SmcqaDfZu7777ru85+m9FAzj9N6kBk04m0n9T/mVptItfu7/1/aY0y633/ctz6NCAatWqmZmrumm3uwZhV6JZW/8/xGbPnm0K/urxdSa5ZuEuNhkCYMICEAIXGxzvb9WqVZ66det6cuXK5SlSpIjn6aef9pw9e/aiA/T9/fDDD562bdt6ChQo4ImJifFUrFjR069fP9/kB53w8MILL5jB2Tlz5vSUKlXKM2bMGN/PP/XUU56rr77aTI544IEHzAB/7ySE06dPex588EFPyZIlTbuKFSvm6d27tyc5Odn38927dzc/r71SOtj6cu677z7zvNmzZ6faf+rUKc8jjzxiXlfPo0ePHp5nnnkm1cDty01YUPv37/e0bt3aTNDQc3zzzTfNOV9uwsLBgwc9LVq0MOeu7dIB6JeasOA/KP1iEzXSDjTX3//kyZPNRBT9vRcqVMgMel+9evVlf0dwN52coO9jIFDZ9D/EsAAAhIfWbdPsbPv27c3KDVda+g1Ii25TAADCSOslauCmY0m7d+/O7x4BI/MGAABgETJvAAAAFiF4AwAAsAjBGwAAgEUI3gAAACxC8AYAAGARgjcAAACLELwBAABYhOANAADAIgRvAAAAFiF4AwAAEHv8HxtIqOSV5VpfAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Create a discrete colormap for binary data\n", "cmap = mcolors.ListedColormap([\"#ffffff\", \"#2ca02c\"]) # Red=False, Green=True\n", "bounds = [0, 0.5, 1]\n", "norm = mcolors.BoundaryNorm(bounds, cmap.N)\n", "\n", "fig, ax = plt.subplots()\n", "im = contains_data.plot(\n", " ax=ax,\n", " cmap=cmap,\n", " norm=norm,\n", " add_colorbar=True,\n", " cbar_kwargs={\n", " \"label\": \"Contains data?\",\n", " \"ticks\": [0.25, 0.75],\n", " },\n", ")\n", "\n", "cbar = im.colorbar\n", "cbar.ax.set_yticklabels([\"False\", \"True\"])\n", "\n", "plt.xlabel(\"Forecast valid time\")\n", "plt.ylabel(\"Forecast lead time [days]\")\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "a005adb2", "metadata": {}, "source": [ "The reason is that forecasts aren't initialized on all days, instead the dataset has gaps between initializations.\n", "\n", "The strategy we're implementing instead (in the `merge_restructure_ds` rule) instead uses `valid_time` as one of the dimensions and introduces a new dimension `sample` that we populate with all forecasts for the same `valid_time`, i.e. all cells that do contain data along each vertical line in the plot above. We compute a look-up table to extract `lead_time` and `init_time` from `valid_time` and `sample`." ] }, { "cell_type": "code", "execution_count": 21, "id": "f773188a", "metadata": { "tags": [ "remove-cell" ] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset> Size: 16GB\n",
       "Dimensions:              (valid_time: 7430, sample: 27, ensemble_member: 11,\n",
       "                          lag: 1, latitude: 18, longitude: 49)\n",
       "Coordinates:\n",
       "  * valid_time           (valid_time) datetime64[ns] 59kB 2003-06-29 ... 2023...\n",
       "  * sample               (sample) int64 216B 0 1 2 3 4 5 6 ... 21 22 23 24 25 26\n",
       "  * ensemble_member      (ensemble_member) int64 88B 0 1 2 3 4 5 6 7 8 9 10\n",
       "  * lag                  (lag) timedelta64[ns] 8B 00:00:00\n",
       "  * latitude             (latitude) float64 144B 30.0 32.5 35.0 ... 70.0 72.5\n",
       "  * longitude            (longitude) float64 392B -80.0 -77.5 ... 37.5 40.0\n",
       "    pressure             int64 8B ...\n",
       "Data variables:\n",
       "    geopotential_height  (valid_time, sample, ensemble_member, lag, latitude, longitude) float64 16GB ...\n",
       "    init_time            (valid_time, sample) datetime64[ns] 2MB ...
" ], "text/plain": [ " Size: 16GB\n", "Dimensions: (valid_time: 7430, sample: 27, ensemble_member: 11,\n", " lag: 1, latitude: 18, longitude: 49)\n", "Coordinates:\n", " * valid_time (valid_time) datetime64[ns] 59kB 2003-06-29 ... 2023...\n", " * sample (sample) int64 216B 0 1 2 3 4 5 6 ... 21 22 23 24 25 26\n", " * ensemble_member (ensemble_member) int64 88B 0 1 2 3 4 5 6 7 8 9 10\n", " * lag (lag) timedelta64[ns] 8B 00:00:00\n", " * latitude (latitude) float64 144B 30.0 32.5 35.0 ... 70.0 72.5\n", " * longitude (longitude) float64 392B -80.0 -77.5 ... 37.5 40.0\n", " pressure int64 8B ...\n", "Data variables:\n", " geopotential_height (valid_time, sample, ensemble_member, lag, latitude, longitude) float64 16GB ...\n", " init_time (valid_time, sample) datetime64[ns] 2MB ..." ] }, "execution_count": 21, "metadata": {}, "output_type": "execute_result" } ], "source": [ "ds_restructured = xr.open_dataset(path_restructured)\n", "ds_restructured.drop_attrs()" ] }, { "cell_type": "markdown", "id": "f46a64c6", "metadata": {}, "source": [ "Repeating the plot above with the new dataset and dimensions `(valid_time, sample)`, we see that we can represent the data much more compactly." ] }, { "cell_type": "code", "execution_count": 22, "id": "825583a0", "metadata": { "tags": [ "remove-input" ] }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAm8AAAHWCAYAAAAhG26oAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjgsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvwVt1zgAAAAlwSFlzAAAPYQAAD2EBqD+naQAAN7ZJREFUeJzt3QucjHX///HPWqx1WinsyuYsK0IoUTnktEk5dlDdFP1yTpJILImt7h9RoXLnVBQdlCLluIVURKTNLW1ZoUWyaK2s+T8+3/9/5j+zdlljd2a+M6/n43Exc10z11zXzOzMe77HMIfD4RAAAABYoZC/DwAAAAB5R3gDAACwCOENAADAIoQ3AAAAixDeAAAALEJ4AwAAsAjhDQAAwCKENwAAAIsUliB39uxZ2b9/v5QqVUrCwsL8fTgAABQIHXP/+PHjUrFiRSlUyDdlM6dOnZLTp0/n2/6KFi0qxYoVy7f9BaugD28a3GJjY/19GAAA+ERqaqpUqlTJJ8GtVHQpOXPsTL7ts3Tp0hITE2PC58CBA82CEAxvWuLmfDPrmwIAgGCUnp5uCiuc33sFTUvcNLhdPeVqKRR56SV9ZzPOyq5hu/i+zoOgD2/OqlINboQ3AECw83UTIQ1u4ZHhPn3MUEeHBQAAAIsQ3gAAACxCeAMAALAI4Q0AAMAihDcAAACLEN4AAAAsQngDAACwCOENAADAIoQ3AAAAixDeAAAALEJ4AwAAsAjhDQAAwCKENwAAAIsQ3gAAACxCeAMAALCIX8NbYmKiNGnSREqVKiXly5eXzp07y65duzxu07t3bwkLC/NYmjZt6rdjBgAACNnwlpSUJAMHDpRNmzbJypUr5cyZM9KuXTs5efKkx+06dOggBw4ccC3Lly/32zEDAAD4U2F/PviKFSs8rs+ZM8eUwG3ZskVuueUW1/qIiAiJjo72wxECAAAEloBq83bs2DHzf9myZT3Wr1u3zoS6WrVqycMPPyxpaWl+OkIAAIAQLnlz53A4ZNiwYXLTTTdJ3bp1Xevj4+OlR48eUrlyZUlJSZExY8ZI69atTemclshll5mZaRan9PR0n50DAABAyIS3QYMGyfbt22X9+vUe6++++27XZQ11jRs3NkFu2bJl0rVr1xw7QYwfP94nxwwAABCS1aaDBw+WpUuXytq1a6VSpUrnvW1MTIwJb7t3785x+6hRo0z1q3NJTU0toKMGAAAIsZI3rSrV4LZkyRLTrq1q1aoXvM+RI0dMINMQlxOtSs2pOhUAACAY+LXkTYcJeeutt2ThwoVmrLeDBw+aJSMjw2w/ceKEDB8+XL766iv59ddfTcDr1KmTXHHFFdKlSxd/HjoAAEDolbzNnDnT/N+yZctzhgzRwXnDw8Nlx44dMn/+fPnrr79MaVurVq1k0aJFJuwBAACEGr9Xm55PZGSkfPbZZz47HgAAgEAXEB0WAAAAkDeENwAAAIsQ3gAAACxCeAMAALAI4Q0AAMAihDcAAACLEN4AAAAsQngDAACwCOENAADAIoQ3AAAAixDeAAAALEJ4AwAAsAjhDQAAwCKENwAAAIsQ3gAAACxCeAMAALAI4Q0AAMAihDcAAACLEN4AAAAsQngDAACwCOENAADAIoQ3AAAAixDeAAAALEJ4AwAAsAjhDQAAwCKENwAAAIsQ3gAAACxCeAMAALAI4Q0AAMAihDcAAACLEN4AAAAsQngDAACwCOENAADAIoQ3AAAAixDeAAAALEJ4AwAAsAjhDQAAwCKENwAAAIsQ3gAAACxCeAMAALAI4Q0AAMAihDcAAACLEN4AAAAsQngDAACwCOENAADAIoQ3AAAAixDeAAAALEJ4AwAAsAjhDQAAwCKENwAAAIsQ3gAAACxCeAMAALAI4Q0AAMAihDcAAACLEN4AAAAsQngDAACwCOENAADAIoQ3AAAAixDeAAAALEJ4AwAAsAjhDQAAwCKENwAAAIsQ3gAAACxCeAMAALCIX8NbYmKiNGnSREqVKiXly5eXzp07y65duzxu43A4ZNy4cVKxYkWJjIyUli1bys6dO/12zAAAACEb3pKSkmTgwIGyadMmWblypZw5c0batWsnJ0+edN3mhRdekClTpsgrr7wi3377rURHR0vbtm3l+PHj/jx0AAAAvyjsn4f9v1asWOFxfc6cOaYEbsuWLXLLLbeYUrepU6fK6NGjpWvXruY28+bNkwoVKsjChQvlkUce8dORAwAA+EdAtXk7duyY+b9s2bLm/5SUFDl48KApjXOKiIiQFi1ayMaNG/12nAAAABLq4U1L2YYNGyY33XST1K1b16zT4Ka0pM2dXnduyy4zM1PS09M9FgAAEDzCwsLOu/Tu3VuCmV+rTd0NGjRItm/fLuvXrz9nm74Q2YNe9nXunSDGjx9fYMcJAAD868CBA67LixYtkrFjx3p0eNQOju7++ecfKVKkiASLgCh5Gzx4sCxdulTWrl0rlSpVcq3XzgkqeylbWlraOaVxTqNGjTLVr84lNTW1gI8eAAD4UnR0tGuJiooyBTrO66dOnZIyZcrI4sWLzQgVxYoVk7feesuMXNGgQQOP/Wi7+ipVqpzT/j4uLs7cr3bt2jJjxoyAe3H9Gt60BE1L3D744ANZs2aNVK1a1WO7XtcXQnuiOp0+fdr0Um3WrFmO+9Q2caVLl/ZYAABAaHnyySdlyJAhkpycLO3bt8/TfWbNmmU6SU6cONHcb9KkSTJmzBjTWTKQ+LXaVIcJ0V6jH330kRnrzVnCpilaizw1SQ8dOtQ8eTVr1jSLXi5evLj07NnTn4cOAAAKQPa26looo8vFGjp0qGukiryaMGGCTJ482XU/LUT68ccf5bXXXpNevXpJoPBreJs5c6b5X4s1sxdZOhsbjhgxQjIyMmTAgAFy9OhRueGGG+Tzzz83YQ8AAASX2NhYj+sJCQmmyvNiNW7c+KJuf+jQIdPUqk+fPvLwww+71usYtFqoFEgK+7va9EK09E1fNG9eOAAAYBcNUO5NnrwpdVMlSpQQd4UKFTond2hHBqezZ8+6qk61oMhdeHi4BJKA6W0KAABQUO3Vy5UrZ5pnuY9YsW3bNtd27Qh55ZVXyi+//CL33XdfQL8QhDcAABD0WrZsaapGddrN7t27m1mePv30U4+gqLV82slB18XHx5uxYzdv3myabelYtIEiIIYKAQAAKEhxcXFm2I/p06dL/fr15ZtvvpHhw4d73KZv377yn//8R+bOnSv16tUzMzrp5eyjYfhbmCMvDc8s77WiDQ11zDeGDQEABCtff985Hy9uZpyER156m7CsjCxJ7p/M93UeUPIGAABgEcIbAACARQhvAAAAFiG8AQAAWITwBgAAYBHCGwAAgEUIbwAAABYhvAEAAFiE8AYAAGARwhsAAIBFCG8AAAAWIbwBAABYhPAGAABgEcIbAACARQhvAAAAFiG8AQAAWITwBgAAYBHCGwAAgEUIbwAAABYhvAEAAFiE8AYAAGARwhsAAIBFCG8AAAAWIbwBAABYhPAGAABgEcIbAACARQhvAAAAFiG8AQAAWITwBgAAYBHCGwAAgEUIbwAAABYhvAEAAFiE8AYAAGARwhsAAIBFCG8AAAAWIbwBAABYhPAGAABgEcIbAACARQhvAAAAFiG8AQAAWITwBgAAYBHCGwAAgEUIbwAAABYhvAEAAFiE8AYAAJDPMjIyZP/+/VIQCG8AAAD5aNWqVVK+fHmJjY2V+Ph4OXHiRH7unvAGAACQn4YPHy7333+/bN68Wc6ePWsCXHp6er7tn5I3AACAfLR7924ZMWKENGzYUD7++GMpXry4XHnllVKtWjXZuXOntGnTxlz2VuH8PFgAAIBQd+WVV5r2blWrVpWiRYvK8uXLZdmyZfLrr7+a6tQuXbrI4cOHvd4/4Q0AACAf9ejRQ1577TVp3ry5uR4eHi533HGHa/vAgQMvaf+ENwAAgHw0duxYSUtLk4JCmzcAAIB8FBERYXqaFhRK3gAAAArIvn37ZOnSpbJ37145ffq0x7YpU6Z4tU/CGwAAQAFYvXq1aeumHRd27doldevWNZ0WHA6HXHfddV7vl2pTAACAAjBq1Ch5/PHH5YcffpBixYrJ+++/L6mpqdKiRQvTqcFbhDcAAIACkJycLL169TKXCxcubKbMKlmypDzzzDPy/PPPe71fwhsAAEABKFGihGRmZprLFStWlD179ri2Mc4bAABAgGnatKls2LBB6tSpIx07djRVqDt27JAPPvjAbPMWHRYAAAAKgPYmdU5KP27cOHN50aJFUqNGDXnxxRe93i/hDQAAoAC4z1+q85vOmDEjX/ZLmzcAAIACCm9Hjhw5Z/1ff/11SRPT+zW8ffHFF9KpUyfTiC8sLEw+/PBDj+29e/c2692XS6kjBgAA8BUd0y0rK+uc9dqJ4ffff7ez2vTkyZNSv359efDBB6Vbt2453qZDhw4yZ84c1/WiRYv68AgBAAAujs6o4PTZZ59JVFSU67qGOR28t0qVKmJleIuPjzfLheYHi46O9tkxAQAAXIrOnTub/7XG0DnOm1ORIkVMcJs8eXLwdlhYt26dlC9fXsqUKWNGJJ44caK5nhstinSOqaLS09N9dKQAAAAiZ8+eNU+DTov17bffyhVXXJGvT0tAd1jQUrkFCxbImjVrTELVJ6B169Ye4Sy7xMREUzzpXGJjY316zAAAAColJSXfg1vAl7zdfffdrss6mWvjxo2lcuXKsmzZMunatWuu84gNGzbMo+SNAAcAAPzVvj8pKUn27t0rp0+f9tg2ZMiQ4Atv2cXExJjwtnv37vO2kdMFAADAn7Zu3Sq33Xab/P333ybElS1b1kyLpWO+aRMwb8NbQFebZqdjpaSmppoQBwAAEMgee+wxMyTan3/+KZGRkbJp0yb57bffpFGjRvK///u/Xu/Xr+FNp4nYtm2bWZx1w3pZixZ12/Dhw+Wrr74y46RoxwV9ArTuuEuXLv48bAAAgAvSTKPzmYaHh5tF2+xrU64XXnhBnnrqKbEyvG3evFkaNmxoFqVt1fTy2LFjzUnq5K133nmn1KpVy3S11f81zJUqVcqfhw0AAHBBOiyIDheiKlSoYAqnlHaodF72aZu3M2fOmNKwPXv2SM+ePU2g2r9/v5QuXVpKliyZp320bNlSHA5Hrtt1YDsAAAAbNWzY0BRUaeFTq1atTOGUtnl78803pV69er4tedP6Wn1QLRUbOHCgHDp0yKzXYkCt6gQAAAh1kyZNcrXTnzBhglx++eXSv39/SUtLk9dff923JW+PPvqoGbbj+++/NwfipG3R+vbt6/XBAAAABIvGjRu7LpcrV06WL1+eL/v1KrytX79eNmzYcM48ozqMx6VMtAoAAIACCG867YNOrJrdvn376EwAAABCup1b2P/rpHAh3333ne/avLVt21amTp3quq4HqUN7JCQkmMHoAAAAQnVS+jvvvNMs7du3Nx07dfIA7aSpS7Fixcw63ebTkrcXX3zR9JqoU6eOnDp1yvQ21VkPdAy2t99+2+uDAQAAsFlCQoLrsvYD0FkUtLNC9tvopAM+DW8VK1Y0A89pUNMiP61G7dOnj9x3331mBGEAAIBQ9+6775qhQrK7//77TWeG2bNn+3acNw1pDz30kFkAAABwblbSTp41a9b0WK/rtPrUW3kOb0uXLs3zTu+44w5vjwcAACAoDB061IzrtmXLFmnatKlZp/ObaombDthb4OFNG+DlhXZeyKknKgAAQCgZOXKkVKtWTaZNmyYLFy406+Li4mTu3Lly1113FXx403ZtAAAAyDsNaZcS1AJuYnoAAAD4KLytXr1abr/9dqlevbrUqFHDXF61apW3uwMAAEBBhbdXXnlFOnToYGZT0HlOdQyT0qVLmwF6dRsAAAAKhldDhSQmJpqBegcNGuRapwGuefPmMnHiRI/1AAAA8HPJW3p6uil5y65du3ZmGwAAADzpaBw6ycHRo0fF5yVvOo7bkiVL5IknnvBY/9FHH0mnTp0u6YBCXb159by6345eO/L9WOD966F4TfIfr0fg4TMLOP84b/Xq1TOzUGlwa9GihWzcuFGKFy8un3zyiZnr1GfhTcco0erRdevWyY033ugadG7Dhg3y+OOPy0svveRRnQoAABBq3nvvPTMVlvr4448lJSVFfvrpJ5k/f76MHj3a5Cafhbc33nhDLrvsMvnxxx/N4lSmTBmzzX3AXsIbAAAIRYcPH5bo6Ghzefny5dKjRw+pVauWKYlzL+jySXjT5AgAAIDcVahQwRRyxcTEyIoVK2TGjBlm/d9//y3h4eHiLa8npgcAAEDuHnzwQTO7goY3rY1s27atWf/1119L7dq1xafhzeFwmHrctWvXSlpa2jlTZ33wwQdeHxAAAEAwGDdunNStW1dSU1NNlWlERIRZr6VuOu+pT8ObDsz7+uuvS6tWrUyRoKZJAAAAeOrevXu2NSK9evWSS+FVeHvrrbdM6ZrOqAAAAELXpp6bzCxLl0rHiY3qHyVNmjQxJVMDBw40i+1Wr15tlpxqKmfPnu278BYVFSXVqlXz6gEBAABy8+233+ZLGAwE48ePl2eeeUYaN27saveWHwp7W4erB6SJMTIyMl8OBAAAIJi8+uqrMnfuXHnggQfydb9ehTdtdPf2229L+fLlpUqVKlKkSBGP7d99911+HR8AAICVTp8+Lc2aNcv3/XoV3nr37i1btmwxowbTYQEAAOBcffv2lYULF8qYMWPE7+Ft2bJl8tlnn8lNN92UrwcDAAAQLE6dOmVG51i1apVce+2159RUTpkyxXfhLTY2NmgaEwIAABSE7du3S4MGDczlH374wWPbpXRe8Cq8TZ48WUaMGGEa4mmbNwAAAHjSyQwKQphDp0u4SDopvc7LdebMGSlevPg5xYB//vmnBAozbkxUlMTNjJPwSO/nEQMAIJBlZWRJcv9kOXbsmE9qx5zfr/n1ePm9v2DmVcnb1KlT8/9IAAAALNe1a1czPIgGUL18Pt5OJ+pVeLvUaR0AAACCUVRUlKs9m14uCF6FN3cZGRnyzz//eKyjuBMAAISiOXPm5Hg5PxXy5k4nT56UQYMGmUF6S5YsadrAuS8AAACQwCl5056m2oNixowZ8q9//UumT58uv//+u7z22mvy3HPP5f9RAgAAWOi9996TxYsXy969e82MC/kxI5VXJW8ff/yxCW7du3eXwoULy8033yxPP/20TJo0SRYsWODVgQAAAASTl156SR588EFTU7l161a5/vrr5fLLL5dffvlF4uPjvd6vV+FNhwKpWrWqq32bc2gQnXHhiy++8PpgAAAAgsWMGTPMDAuvvPKKFC1a1NRcrly5UoYMGWKGRPFpeKtWrZr8+uuv5nKdOnVMcaCzRK5MmTJeHwwAAECw2Lt3r2ti+sjISDl+/Li5/MADD8jbb7/t2/CmRYDff/+9uTxq1CiTLCMiIuSxxx6TJ554wuuDAQAACBbR0dFy5MgRc7ly5cqyadMmczklJUW8mCPh0josaEhzatWqlfz000+yefNmqV69utSvX9/rgwEAAAgWrVu3NrWS1113nfTp08fkJ+3AoJnpQgP45lt4+/rrr037NvdGdvPnz5eEhAQzfEjnzp3l5ZdfNqVwAAAAoez111+Xs2fPmsv9+vWTsmXLyvr166VTp07muk+qTceNGyfbt293Xd+xY4dJkm3atDHVp5ouExMTvT4YAACAYLFv3z4JD///86rfddddpgfq4MGD5eDBg74Jb9u2bZNbb73Vdf2dd96RG264QWbNmmWKAvWAnJ0XAAAAQlnVqlXl0KFD5x21o8DD29GjR6VChQqu60lJSdKhQwfX9SZNmkhqaqrXBwMAABAsHA6Ha55TdydOnJBixYr5ps2bBjftIREbG2tGCdaRgcePH+/arl1gixQp4vXBAAAA2G7YsGHmfw1uY8aMkeLFi7u2ZWVlmT4EDRo08E1401K2kSNHyvPPPy8ffvihORidXcFJ28Npj1MAAIBQtXXrVlfJm/YP0AF6nfSyjswxfPhw34S3Z5991nRtbdGihZmQft68eR4HNHv2bGnXrp3XBwMAAGC7tWvXusbFnTZtmpmNKj9dVHgrV66cfPnll2ZKBw1v7j0o1LvvvmvWAwAAhLo5c+YUyH69GqQ3Kioqx/U6fgkAAADEjIH73HPPyerVqyUtLc015puTTlDvs/AGAACA8+vbt68ZmUPnMo2Jicmx56k3CG8AAAAF4NNPP5Vly5ZJ8+bN83W/Xk1MDwAAgPO77LLLCqRJGeENAACgAEyYMEHGjh0rf//9d77ul2pTAACAAjB58mTZs2ePmeSgSpUq50xkoJMdeIPwBgAAUAA6d+5cELslvAEAABSEhISEAtkvJW8AAAAFaMuWLZKcnGyGCqlTp440bNjwkvZHeAMAACgAOjDvPffcI+vWrZMyZcqYuU51lqpWrVrJO++8Y2au8ga9TQEAAArA4MGDJT09XXbu3Cl//vmnHD16VH744QezbsiQIV7vl5I3AACAArBixQpZtWqVxMXFudZpten06dOlXbt2Xu+XkjcAAIACoHOZZh8eROm67POcXgzCGwAAQAFo3bq1PProo7J//37Xut9//10ee+wxufXWW+0Mb1988YV06tRJKlasaHpgfPjhhx7btWHfuHHjzPbIyEhp2bKlqTcGAAAIdK+88oocP37cDNBbvXp1qVGjhlStWtWse/nll+1s83by5EmpX7++PPjgg9KtW7dztr/wwgsyZcoUmTt3rtSqVUueffZZadu2rezatUtKlSrll2MGAADIi9jYWDOLwsqVK+Wnn34yhVLa5q1NmzZyKfwa3uLj482SEz3BqVOnyujRo6Vr165m3bx588wUEwsXLpRHHnnEx0cLAABw8bTgSZf8ErBt3lJSUuTgwYMevTEiIiKkRYsWsnHjRr8eGwAAQG7WrFljSth0SJDsdJy3a665Rr788ksJuvCmwU1pSZs7ve7clpPMzEzzZLkvAAAAvqI1hw8//LCULl36nG1RUVGm9lCbhQVdeHPSjgzZq1Ozr3OXmJhonhjnovXNAAAAvvL9999Lhw4dct2utYo6ZVbQhbfo6Gjzf/ZSNp1qIntpnLtRo0aZIknnkpqaWuDHCgAA4PTHH3/kOL6bU+HCheXQoUMSdOFNu9JqgNMeGk6nT5+WpKQkadasWa7303ZxWkzpvgAAAPjKlVdeKTt27Mh1+/bt2yUmJsbO8HbixAnZtm2bWZydFPTy3r17TdXo0KFDZdKkSbJkyRIzF1jv3r2lePHi0rNnT38eNgAAQK5uu+02GTt2rJw6deqcbRkZGZKQkCC33367WDlUyObNm6VVq1au68OGDTP/9+rVy4ztNmLECHOSAwYMMJO53nDDDfL5558zxhsAAAhYTz/9tHzwwQdmjNpBgwbJ1VdfbQqlkpOTzbymWVlZZig0b4U5tAdAENPeptpxIW5mnIRHhvv7cAAAKBBZGVmS3D/ZtPf2RZMh5/drfj1efu/P33777Tfp37+/fPbZZ6azpdIA1759e5kxY4aZdcHKkjcAAIBgVLlyZVm+fLmpOfz5559NgKtZs6Zcdtlll7xvwhsAAEAB0bDWpEmTfN1nwPY2BQAAwLkIbwAAABYhvAEAAFiE8AYAAGARwhsAAIBFCG8AAAAWIbwBAABYhPAGAABgEcIbAACARQhvAAAAFiG8AQAAWITwBgAAYBHCGwAAgEUIbwAAABYhvAEAAFiE8AYAAGARwhsAAIBFCG8AAAAWIbwBAABYhPAGAABgEcIbAACARQhvAAAAFiG8AQAAWITwBgAAYBHCGwAAgEUIbwAAABYhvAEAAFiE8AYAAGARwhsAAIBFCG8AAAAWIbwBAABYhPAGAABgEcIbAACARQhvAAAAFiG8AQAAWITwBgAAYBHCGwAAgEUIbwAAABYhvAEAAFiE8AYAAGARwhsAAIBFCG8AAAAWIbwBAABYhPAGAABCxty5c6VMmTJiM8IbAACwTu/evSUsLOyc5eeff5ZgV9jfBwAAAOCNDh06yJw5czzWlStXLuifTEreAACAlSIiIiQ6OtpjmTZtmtSrV09KlCghsbGxMmDAADlx4kSu+/j++++lVatWUqpUKSldurQ0atRINm/e7Nq+ceNGueWWWyQyMtLsb8iQIXLy5EnxJ8IbAAAIGOnp6R5LZmbmRd2/UKFC8tJLL8kPP/wg8+bNkzVr1siIESNyvf19990nlSpVkm+//Va2bNkiI0eOlCJFiphtO3bskPbt20vXrl1l+/btsmjRIlm/fr0MGjRI/IlqUwAAEDC0dMtdQkKCjBs3LsfbfvLJJ1KyZEnX9fj4eHn33Xdd16tWrSoTJkyQ/v37y4wZM3Lcx969e+WJJ56Q2rVrm+s1a9Z0bfv3v/8tPXv2lKFDh7q2aTBs0aKFzJw5U4oVKyb+QHgDAAABIzU11VRfuleN5qZVq1YmRDlpVenatWtl0qRJ8uOPP5qSuzNnzsipU6dMVaduz27YsGHSt29fefPNN6VNmzbSo0cPqV69utmmJXHaAWLBggWu2zscDjl79qykpKRIXFyc+APVpgAAIGBocHNfzhfeSpQoITVq1HAtp0+flttuu03q1q0r77//vglf06dPN7f9559/ctyHlurt3LlTOnbsaKpY69SpI0uWLDHbNKQ98sgjsm3bNteibeR2797tCnj+QMkbAAAICps3bzYlbZMnTzZt39TixYsveL9atWqZ5bHHHpN7773X9GDt0qWLXHfddSbYaTAMJJS8AQCAoFC9enUT3l5++WX55ZdfTFXoq6++muvtMzIyTOeDdevWyW+//SYbNmwwHRec1aFPPvmkfPXVVzJw4EBT6qYlbkuXLpXBgweLPxHeAABAUGjQoIFMmTJFnn/+eVN1qm3VEhMTc719eHi4HDlyRP71r3+Zkre77rrLdHoYP3682X7ttddKUlKSCW0333yzNGzYUMaMGSMxMTHiT2EObXkXxLSxYlRUlMTNjJPwyHB/Hw4AAAUiKyNLkvsny7Fjxzwa/Bf092t+PV5+7y+YUfIGAABgEcIbAACARQhvAAAAFiG8AQAAWITwBgAAYBHCGwAAgEUIbwAAABYJ6PCm842FhYV5LNHR0f4+LAAAAL8J+LlNr7nmGlm1apXHaMgAAAChKuDDW+HChSltAwAAsKHaVOl8YhUrVpSqVavKPffcYyaaBQAACFUBXfJ2ww03yPz5881ksX/88Yc8++yz0qxZM9m5c6dcfvnlOd4nMzPTLO5zpQEAAASLgC55i4+Pl27dukm9evWkTZs2smzZMrN+3rx5ud4nMTHRTGzrXGJjY314xAAAACEc3rIrUaKECXJalZqbUaNGybFjx1xLamqqT48RAAAgZKtNs9Pq0OTkZLn55ptzvU1ERIRZAAAAglFAl7wNHz5ckpKSJCUlRb7++mvp3r27acPWq1cvfx8aAACAXwR0ydu+ffvk3nvvlcOHD0u5cuWkadOmsmnTJqlcubK/Dw0AAMAvAjq8vfPOO/4+BAAAgIAS0NWmAAAA8ER4AwAAsAjhDQAAwCKENwAAAIsQ3gAAACxCeAMAALAI4Q0AAMAihDcAAACLEN4AAAAsQngDAACwCOENAADAIoQ3AAAAixDeAAAALEJ4AwAAsAjhDQAAwCKENwAAAIsQ3gAAACxCeAMAALAI4Q0AAMAihDcAAACLEN4AAAAsQngDAACwCOENAADAIoQ3AAAAixDeAAAALEJ4AwAAsAjhDQAAwCKENwAAAIsQ3gAAACxCeAMAALAI4Q0AAMAihDcAAACLEN4AAAAsQngDAACwCOENAADAIoQ3AAAAixDeAAAALEJ4AwAAsAjhDQAAwCKENwAAAIsQ3gAAACxCeAMAALAI4Q0AAMAihDcAAACLEN4AAAAsQngDAACwCOENAADAIoQ3AAAAixDeAAAALEJ4AwAAsAjhDQAAwCKENwAAAIsQ3gAAACxCeAMAALAI4Q0AAMAihDcAAACLEN4AAAAsQngDAACwCOENAADAIoQ3AAAAixDeAAAALEJ4AwAAsAjhDQAAwCKENwAAAIsQ3gAAACxiRXibMWOGVK1aVYoVKyaNGjWSL7/80t+HBAAA4BcBH94WLVokQ4cOldGjR8vWrVvl5ptvlvj4eNm7d6+/Dw0AAMDnAj68TZkyRfr06SN9+/aVuLg4mTp1qsTGxsrMmTP9fWgAAAA+F9Dh7fTp07JlyxZp166dx3q9vnHjxhzvk5mZKenp6R4LAABAsAjo8Hb48GHJysqSChUqeKzX6wcPHszxPomJiRIVFeVatJQOAAAgWAR0eHMKCwvzuO5wOM5Z5zRq1Cg5duyYa0lNTfXRUQIAABS8whLArrjiCgkPDz+nlC0tLe2c0jiniIgIswAAAASjgC55K1q0qBkaZOXKlR7r9XqzZs38dlwAAAD+EtAlb2rYsGHywAMPSOPGjeXGG2+U119/3QwT0q9fP38fGgAAgM8FfHi7++675ciRI/LMM8/IgQMHpG7durJ8+XKpXLmyvw8NAADA5wI+vKkBAwaYBQAAINRZEd4uhfZMVVkZWf4+FAAACozze875vecr+TWeKuOy5l2Yw9evso/t27ePsd4AACFDh8iqVKlSgT/OqVOnzLzjuY276o3SpUtLTEyMFCpUSAYOHGgWhGB4O3v2rOzfv19KlSp1zthwmvJ1EF99o+sbJhiFwjmGynlyjsGD1zI4BNrrqF/nx48fl4oVK5rw4wsa4HQ2pPwcZaJYsWL5tr9gFfTVpvoGvtAvEP2jC4Q/vIIUCucYKufJOQYPXsvgEEivo84s5EsatAhbvhfQ47wBAADAE+ENAADAIiEd3nQarYSEhKCeTisUzjFUzpNzDB68lsEhFF5HBKag77AAAAAQTEK65A0AAMA2hDcAAACLEN4AAAAsEtThrXfv3tK5c2cJFTNmzDCjXeuYO40aNZIvv/zS47nQQYrdl6ZNm0qwnecff/xhzlUHqSxevLh06NBBdu/eLTb54osvpFOnTuYc9HX68MMPPbaPGzdOateuLSVKlJDLLrtM2rRpI19//bUE0zlmf686l3//+99ii8TERGnSpIkZILx8+fLms2jXrl0et/nggw+kffv2csUVV5jz27Ztm9gkL+do+2dPXs4xGD53YJegDm+hZNGiRTJ06FAZPXq0bN26VW6++WaJj4+XvXv3um6jHygHDhxwLcuXL5dgOk/te6MfrL/88ot89NFHZnvlypVNuDl58qTYQo+1fv368sorr+S4vVatWmbbjh07ZP369VKlShVp166dHDp0SILlHN3fp7rMnj3bfOl369ZNbJGUlGSm9tm0aZOsXLlSzpw5Y14n9/eiXm7evLk899xzYqO8nKPtnz0XOsdg+dyBZRxBrFevXo4777zTXP70008dzZs3d0RFRTnKli3r6Nixo+Pnn3923TYlJUV73Tref/99R8uWLR2RkZGOa6+91rFx40aHDa6//npHv379PNbVrl3bMXLkyHOeC5ud7zx37dplXsMffvjBte3MmTPm9Z41a5bDRno+S5YsOe9tjh07Zm63atUqR7Ceo753W7du7bBZWlqaOdekpKRztjk/f7Zu3eoItnMMls+e3M4xGD93EPhCpuRNfwENGzZMvv32W1m9erWZNqtLly5m7lN3WqIzfPhwU32hJRz33nuv+aUVyHReuS1btphfg+70+saNG13X161bZ4r99bwefvhhSUtLE5tc6DwzMzPNdfepWsLDw81ceVpCFYz0OXn99dfNlDhakhWMtEpq2bJl0qdPH7HZsWPHzP9ly5aVYJXbOdr+2XO+cwzFzx34X8iEN61u6dq1q9SsWVMaNGggb7zxhql2+vHHHz1up8GtY8eO5kNm/Pjx8ttvv8nPP/8sgezw4cOSlZUlFSpU8Fiv1w8ePGgua9XiggULZM2aNTJ58mQTYlu3bu364LHBhc5T24FpdcWoUaPk6NGjJthodZRu06qaYPLJJ59IyZIlzRfGiy++aKpztN1UMJo3b55pb6R/v7bSAkb98XjTTTdJ3bp1JRjldo7B8NlzvnMMpc8dBI6gn5jeac+ePTJmzBjTbkFDgLPETdtKuX/QXHvtta7LMTEx5n/9lah/oIFO2wRl/6Bxrrv77rtd6/V8GzdubD5wtETDti/F3M6zSJEi8v7775sSGv1VrL9+td2JfnkEm1atWpnSYX0vz5o1S+666y7TaUFLN4KNtne77777rJ78etCgQbJ9+/agLonJ7RyD6bMnp3MMpc8dBI6QKXnTnm1HjhwxX3T6Jefsnae/ktzpH2L2kJC9ajXQaImLfmA4S9mcNHRmL6VyD6b6AWpTj6i8nKf2PtVQ89dff5lfvStWrDCvu/ZODSba07RGjRqm156WIhcuXNj8H2y0J7H27Ovbt6/YavDgwbJ06VJZu3atVKpUSYLRxZyjjZ89FzrHUPncQeAIifCmf0TJycny9NNPy6233ipxcXGmeDtYaNsK/fDQqjN3er1Zs2a5Piepqamu0sVgO09tA1auXDnzBbF582a58847JZhp6aON1VAXooFUX3Mb2/Ppa6IlNTociFYZBuMXuTfnaNtnz8WcY6h97sB/QqLaVMfCuvzyy03Dbv3A0KrSkSNHSjDRdhgPPPCAqZK48cYbzbnqefbr109OnDhhxgbTdn96/r/++qs89dRTpiRLO20Ey3mqd99913x4XnXVVaZN46OPPmq68Wfv5BDI9PVyb2eZkpJiftVrlYy+jydOnCh33HGHeS31i1DHvdu3b5/06NFDguEc9bVT6enp5vXUdlI20uElFi5caIaP0DZ7zhJj/YKPjIw0l//880/z/t2/f7+57hw/LDo62iy2n2MwfPbk5XUMhs8dWMYRxB544AFHt27dzOWVK1c64uLiHBEREWYIkHXr1nkMUZBTV/2jR4+adWvXrnXYYPr06Y7KlSs7ihYt6rjuuutcXdn//vtvR7t27RzlypVzFClSxHHVVVeZ7vt79+512Ci381TTpk1zVKpUyXWeTz/9tCMzM9NhE32/6fsu+6KvWUZGhqNLly6OihUrmvOPiYlx3HHHHY5vvvnGESzn6PTaa6+ZIXv++usvh41yOj9d5syZ47qNXs7pNgkJCY5gOMdg+OzJy+sYDJ87sEuY/iNBSgeG1HZBuQ0ECgAAYJugbPOm7dm0J5OOLaS9fgAAAIJFULZ5e+ihh8xYQo8//jgNRgEAQFAJ6mpTAACAYBOU1aYAAADBivAGAABgEcIbAACARQhvAAAAFiG8AQAAWITwBiDg9O7d20wv5NSyZUsZOnToee9TpUoVmTp16kU9jk7XFBYWZqbmQnBJTEyUJk2amCmtypcvb95PzunHnHSwBZ2+q2LFimaqK32f7dy50+M2OgWfri9durR5r+jk8zmNLarT9umUWbro5Zxul50+doMGDfLhbBFqCG9AAYUP/aDPvrjP52mbuXPnSpkyZfzy2Dop+IQJE/I1EKrY2Fg5cOCA1K1b9xKPEIEmKSnJzEu6adMmWblypZw5c8bMNXry5EnXbV544QWZMmWKmYVHxwbV+WTbtm0rx48fd93m77//NrP16JysuenZs6f5AbBixQqz6GUNcEBBIbwBBUQ/8DUYuC9Vq1b1al+nT5+WUKYT1msJSn4LDw83X9iFCwfleOUhTUOUBvZrrrlG6tevL3PmzJG9e/fKli1bXKVuWlI7evRo6dq1qwnw8+bNM2FNJ6J30hLfkSNHStOmTXN8nOTkZPNY//nPf+TGG280y6xZs+STTz45p6TvQjRAani84oorTAleixYt5LvvvvO4jf4I1Mfq0qWLFC9eXGrWrClLly716jmCvQhvQAGJiIgwwcB90bDgLBW4/vrrzW1iYmLMl4OWDDhpNc2gQYNk2LBh5oNcP9DVjz/+KLfddpuULFlSKlSoYH7dHz582HW/s2fPyvPPP2/m9NV9X3XVVTJx4kTX9ieffFJq1aplPvSrVasmY8aMkX/++ce1/fvvv5dWrVqZoKTVRI0aNZLNmzebqeYefPBBOXbsmKsUUat8stMvK932008/eazX0g2t1tQvzKysLOnTp48JslpVdfXVV8u0adPO+1xmrzZNS0uTTp06mfvrfhYsWHDe++ux6hfzRx995Dp+Pafs1aa6Tq9/9tln0rBhQ7P/1q1bm8f79NNPJS4uzjwv9957r/mSd9Lz0lIcfU71PhoW3nvvvfMeE3xL37vOHwIqJSVFDh48aErjnPRvRgPTxo0b87zfr776ygStG264wbVOg56uu5j9KC3x69Wrl3z55ZemxFCDmf69u5cEqvHjx8tdd90l27dvN9vvu+8++fPPPy/qsWA3fm4CPvb777+bD1wtFZg/f74JOg8//LAUK1bMIxBp2Ojfv79s2LDBhAMtudMvFr2thqGMjAwTxvRDfM2aNeY+o0aNMr/6X3zxRbnpppvMfdyDlIYyrf7UNj47duww+9J1I0aMMNv1S0BDy8yZM03Q1FBTpEgRadasmSmlGDt2rKs0QQNkdhrENPBpmHKv5tSSDK1a0mCkAbNSpUqyePFiE0z1C+5//ud/TIjVc8kLfe5SU1PNeRctWlSGDBliAlZuhg8fbkpI0tPTTQmM80t8//79Od5eXwetStOQq8eki36x63mcOHHClHq8/PLL5vlXTz/9tKna1edNv3C/+OILuf/++6VcuXLmNYN/6d+P/hDSvwlnFbkGN6U/gtzp9d9++y3P+9b9aJu67HSd8zHySn8ouHvttdfksssuMz/2br/9do/3v/6AUJMmTTLvxW+++caU9iNE6PRYAPJXr169HOHh4Y4SJUq4lu7du5ttTz31lOPqq692nD171nX76dOnO0qWLOnIysoy11u0aOFo0KCBxz7HjBnjaNeunce61NRUnd7OsWvXLkd6erojIiLCMWvWrDwf5wsvvOBo1KiR63qpUqUcc+fOzfG2c+bMcURFRV1wn1OmTHFUq1bNdV2PTY9x586dud5nwIABjm7dunk8f3feeafruj4fjz76qMf+Nm3a5NqenJxs1r344ou5Pkb2faqUlBRzv61bt5rra9euNddXrVrluk1iYqJZt2fPHte6Rx55xNG+fXtz+cSJE45ixYo5Nm7c6LHvPn36OO69995cjwe+o++vypUrm78Xpw0bNpjXdf/+/R637du3r+u1ded8bxw9etRj/cSJEx21atU65/Y1atQw7x1Vp04d1+dAhw4dXLdJSEhw1K9f33X9jz/+MO+tmjVrOkqXLm1uHxYWZj4fnPQYFi9e7PFYett58+Zd5LMCm1HyBhQQrX7UkhinEiVKmP+1BEjbxWgplFPz5s1Nic6+fftMVadq3Lixx/60rc7atWtzLPHas2eP6d2WmZkpt956a67HpFV5WoKmHSf08bSqVqsBnbR0om/fvvLmm29KmzZtpEePHlK9evWLOu977rlHnnjiCVPto9VHWgqnPerq1Knjus2rr75q2u1oCYeWIGqbvrz2utPnT9uouT8/tWvXztfOFNdee61HSYyzmtl9nZZ0OKuyT5065aradtJz0lJM+NfgwYNNmzAtDdUSXydtxqC0dExLfZ20BDd7adz56H7++OOPc9YfOnTItZ/ly5e7midotXputERN76d/o5UrVzalvfpZkb3Nq5aGu3OWaCN00OYNKCAa1rTtmXNxfkHoj2f34OZcp9zXO8Oek344azsvrcp0X3bv3i233HLLeb8UlIYpDVbx8fGmMfXWrVtNY233LwatLtShEjp27GiqJDVwLVmy5KLOW89Tg6uz0ffbb79tqhCdtLr0sccek4ceekg+//xzcw7ani6vnTJyeq7ym/uXoz7O+b4snf8vW7bM43XRUEe7N//R94m2G9XqbH0vZ+8spNc1eGlPVCd9D2oVpTYTyCsNV9qezhnm1ddff23WOfejQcz5OXDllVfmui9t66ZNALRZhXa00PDm3qYVcKLkDfAxDUTvv/++R4jTdl/a9ux8H+zXXXeduZ82/M+pd6S2tdIAt3r1alN6lp22ndMvEQ1sTjm17dEODbpowNJ2NdpGTNt4adsy7WyQF9p2TtuD6f21VFBDo/sXlH6pDRgwwLVOb5NX2mlASwy1I4V2+lDaDu9C42pdzPFf7OupX7Lak5H2bYFDhwnRHxDaSUX/tpztz7Qjgf6d6N+edoLRNmP6t6OLXtZSVm2f6aT308U5zI+2FdX9aQm5tpvU96O2NdP2o9pGTWkbTm2jpm1AL4aGOy311lJlbZ+pJdgX+lGG0ETJG+BjGlq0sb1W52hnAv1ySUhIMFWWhQoVOu+XkfYo00Ckv/J/+eUXU3KlJVgaSrTDgwYm7XygHSE0EGlp2xtvvOH6YtCA8c4775htL730kkepmlZfakmF9rjUUKdhT4cu0C8npaFRq1o1HGppgHtvy+x06AX98tEOF1oK5x5K9Tg0eGmPzv/+97+mx6s+Tl7pF6Lzy1JLOLQ6WcPqhb7k9Pi1d54GPT1+9162l0K/yLVDhIZd7WSiz62Wak6fPt1ch39okwUt/dKeyloa7FwWLVrkuo3+rWiA079JDUzamUj/ptyHpdEqfq3+1veb0lJuve4+PIc2DahXr57puaqLVrtrCLsQLbV1/yE2e/ZsM+Cv7l97kmspXE6dIQA6LAAFIKfG8e7WrVvnaNKkiaNo0aKO6Ohox5NPPun4559/cmyg7+6///2vo0uXLo4yZco4IiMjHbVr13YMHTrU1flBOzw8++yzpnF2kSJFHFdddZVj0qRJrvs/8cQTjssvv9x0jrj77rtNA39nJ4TMzEzHPffc44iNjTXHVbFiRcegQYMcGRkZrvv369fP3F9rpbSx9fn06NHD3G727Nke60+dOuXo3bu3eVw9j/79+ztGjhzp0XD7fB0W1IEDBxwdO3Y0HTT0HOfPn2/O+XwdFtLS0hxt27Y1567HpQ3Qc+uw4N4oPaeOGtkbmuvzP23aNNMRRZ/3cuXKmUbvSUlJ532OENq0c4K+j4GLFab/kGEBAPANHbdNS2e7d+9uZm640NRvQHZUmwIA4EM6XqIGN21L2q9fP557XDRK3gAAACxCyRsAAIBFCG8AAAAWIbwBAABYhPAGAABgEcIbAACARQhvAAAAFiG8AQAAWITwBgAAYBHCGwAAgEUIbwAAAGKP/wOV+E7PQPrI0gAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "# Create a discrete colormap for binary data\n", "cmap = mcolors.ListedColormap([\"#ffffff\", \"#2ca02c\"]) # Red=False, Green=True\n", "bounds = [0, 0.5, 1]\n", "norm = mcolors.BoundaryNorm(bounds, cmap.N)\n", "\n", "fig, ax = plt.subplots()\n", "im = (\n", " (~np.isnat(ds_restructured.init_time.sel(valid_time=range_vt)))\n", " .transpose()\n", " .plot(\n", " ax=ax,\n", " cmap=cmap,\n", " norm=norm,\n", " add_colorbar=True,\n", " cbar_kwargs={\n", " \"label\": \"Contains data?\",\n", " \"ticks\": [0.25, 0.75],\n", " },\n", " )\n", ")\n", "\n", "\n", "cbar = im.colorbar\n", "cbar.ax.set_yticklabels([\"False\", \"True\"])\n", "\n", "plt.title(\"\")\n", "plt.xlabel(\"Forecast valid time\")\n", "plt.ylabel(\"Sample\")\n", "plt.tight_layout()\n", "plt.show()\n" ] }, { "cell_type": "markdown", "id": "8861d9c6", "metadata": {}, "source": [ "The size of the `sample` dimension is chosen as the maximum number of states with same `valid_time`. This is already smaller than the lead_time dimension:" ] }, { "cell_type": "code", "execution_count": 23, "id": "1a9ffdc6", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Length lead_time dimension: 46\n", "Length sample dimension: 27\n" ] } ], "source": [ "print(\"Length lead_time dimension:\", len(ds_combined.lead_time))\n", "print(\"Length sample dimension:\", len(ds_restructured.sample))" ] }, { "cell_type": "markdown", "id": "b5569cbf", "metadata": {}, "source": [ "In our dataset, everything larger than `sample == 13` contains almost exclusively missing data. Therefore, we usually restrict the data to samples `0` to `13` by setting `n_samples: 14` in the corresponding configuration file. This is especially important as the computational demand for precomputing similarities scales quadratically with the number of included states." ] }, { "cell_type": "markdown", "id": "bf2da704", "metadata": {}, "source": [ "### Splitting valid time dimension into year and dayofyear" ] }, { "cell_type": "markdown", "id": "cfdbc01e", "metadata": {}, "source": [ "What we want to compute is similarities between states with a similar *day of year* of `valid_time`, not with a similar `valid_time` itself.\n", "\n", "The next step (rule `to_year_dayofyear_format`) therefore consists of splitting the `valid_time`dimension into two separate dimensions `year` and `dayofyear`. We chunk the dataset along the `dayofyear` dimension to allow memory-efficient computations:" ] }, { "cell_type": "code", "execution_count": 24, "id": "6c038876", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'geopotential_height' (dayofyear: 366, year: 21, sample: 27,\n",
       "                                         ensemble_member: 11, lag: 1,\n",
       "                                         latitude: 18, longitude: 49)> Size: 16GB\n",
       "dask.array<open_dataset-geopotential_height, shape=(366, 21, 27, 11, 1, 18, 49), dtype=float64, chunksize=(1, 21, 27, 11, 1, 18, 49), chunktype=numpy.ndarray>\n",
       "Coordinates:\n",
       "  * dayofyear        (dayofyear) int64 3kB 1 2 3 4 5 6 ... 362 363 364 365 366\n",
       "  * year             (year) int64 168B 2003 2004 2005 2006 ... 2021 2022 2023\n",
       "  * sample           (sample) int64 216B 0 1 2 3 4 5 6 ... 20 21 22 23 24 25 26\n",
       "  * ensemble_member  (ensemble_member) int64 88B 0 1 2 3 4 5 6 7 8 9 10\n",
       "  * lag              (lag) timedelta64[ns] 8B 00:00:00\n",
       "  * latitude         (latitude) float64 144B 30.0 32.5 35.0 ... 67.5 70.0 72.5\n",
       "  * longitude        (longitude) float64 392B -80.0 -77.5 -75.0 ... 37.5 40.0
" ], "text/plain": [ " Size: 16GB\n", "dask.array\n", "Coordinates:\n", " * dayofyear (dayofyear) int64 3kB 1 2 3 4 5 6 ... 362 363 364 365 366\n", " * year (year) int64 168B 2003 2004 2005 2006 ... 2021 2022 2023\n", " * sample (sample) int64 216B 0 1 2 3 4 5 6 ... 20 21 22 23 24 25 26\n", " * ensemble_member (ensemble_member) int64 88B 0 1 2 3 4 5 6 7 8 9 10\n", " * lag (lag) timedelta64[ns] 8B 00:00:00\n", " * latitude (latitude) float64 144B 30.0 32.5 35.0 ... 67.5 70.0 72.5\n", " * longitude (longitude) float64 392B -80.0 -77.5 -75.0 ... 37.5 40.0" ] }, "execution_count": 24, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xr.open_zarr(path_dayofyear_year)[\"geopotential_height\"].drop_attrs()" ] }, { "cell_type": "markdown", "id": "47725df0", "metadata": {}, "source": [ "### A note on the provided ERA5 datasets\n", "\n", "In order to make it easier to use ERA5 data within the *unseen-awg* framework, we cast it to the same format as the reforecast data. This means, that we introduce \"dummy dimensions\" `init_time` and `ensemble_member` with size 1 and convert the ERA5 time axis to the `lead_time` with respect to the first included time step:" ] }, { "cell_type": "code", "execution_count": 25, "id": "697ad0f9", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'geopotential_height' (ensemble_member: 1, init_time: 1,\n",
       "                                         lead_time: 7670, latitude: 18,\n",
       "                                         longitude: 49)> Size: 27MB\n",
       "[6764940 values with dtype=float32]\n",
       "Coordinates:\n",
       "  * ensemble_member  (ensemble_member) int64 8B 0\n",
       "  * init_time        (init_time) datetime64[ns] 8B 2003-01-01\n",
       "  * lead_time        (lead_time) timedelta64[ns] 61kB 0 days ... 7669 days\n",
       "  * latitude         (latitude) float64 144B 30.0 32.5 35.0 ... 67.5 70.0 72.5\n",
       "  * longitude        (longitude) float64 392B -80.0 -77.5 -75.0 ... 37.5 40.0\n",
       "    pressure         int64 8B ...
" ], "text/plain": [ " Size: 27MB\n", "[6764940 values with dtype=float32]\n", "Coordinates:\n", " * ensemble_member (ensemble_member) int64 8B 0\n", " * init_time (init_time) datetime64[ns] 8B 2003-01-01\n", " * lead_time (lead_time) timedelta64[ns] 61kB 0 days ... 7669 days\n", " * latitude (latitude) float64 144B 30.0 32.5 35.0 ... 67.5 70.0 72.5\n", " * longitude (longitude) float64 392B -80.0 -77.5 -75.0 ... 37.5 40.0\n", " pressure int64 8B ..." ] }, "execution_count": 25, "metadata": {}, "output_type": "execute_result" } ], "source": [ "xr.open_dataarray(path_circulation_ds_era5).drop_attrs()" ] }, { "cell_type": "markdown", "id": "e5d7f68a", "metadata": {}, "source": [ "The rest of the processing pipeline follows what is done for reforecast data." ] } ], "metadata": { "kernelspec": { "display_name": "unseen-awg", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.12" } }, "nbformat": 4, "nbformat_minor": 5 }