diff --git a/CHANGELOG.rst b/CHANGELOG.rst index 9979f7af4e0bc5eda18a9c09ae1cb211402e4734..b163242b0757091d00f95816694e1e58704aecc0 100644 --- a/CHANGELOG.rst +++ b/CHANGELOG.rst @@ -17,6 +17,7 @@ CHANGELOG * Documentation * Removed negative values on the horizontal interpolation due to unmapped NaNs values. * Improved load_nes.py removing redundant code + * New functionalities for vertical extrapolation. (`#74 `_) * Bugfix: * Vertical interpolation for descendant level values (Pressure) (`#71 `_) diff --git a/nes/methods/vertical_interpolation.py b/nes/methods/vertical_interpolation.py index 1605ed7ddfc07f14b0b7b5ad40c5cae0a199f45a..7260fdbe191484129649b0feff7e3efcc4d95728 100644 --- a/nes/methods/vertical_interpolation.py +++ b/nes/methods/vertical_interpolation.py @@ -26,7 +26,70 @@ def add_4d_vertical_info(self, info_to_add): return None -def interpolate_vertical(self, new_levels, new_src_vertical=None, kind='linear', extrapolate=None, info=None, +def parse_extrapolate(extrapolate) -> tuple: + """ + Parses the 'extrapolate' parameter and returns a tuple representing the extrapolation options. + + Parameters + ---------- + extrapolate : bool or tuple or None or number or NaN + If bool: + - If True, both extrapolation options are set to 'extrapolate'. + - If False, extrapolation options are set to ('bottom', 'top'). + If tuple: + - The first element represents the extrapolation option for the lower bound. + - The second element represents the extrapolation option for the upper bound. + - If any element is bool: + - If True, it represents 'extrapolate'. + - If False: + - If it's the first element, it represents 'bottom'. + - If it's the second element, it represents 'top'. + - If any element is None, it is replaced with numpy.nan. + - Other numeric values are kept as they are. + - If any element is NaN, it is kept as NaN. + If None: + - Both extrapolation options are set to (NaN, NaN). + If number: + - Both extrapolation options are set to the provided number. + If NaN: + - Both extrapolation options are set to NaN. + + Returns + ------- + tuple + A tuple representing the extrapolation options. If the input is invalid, it returns + ('extrapolate', 'extrapolate'). + """ + if isinstance(extrapolate, bool): + if extrapolate: + extrapolate_options = ('extrapolate', 'extrapolate') + else: + extrapolate_options = ('bottom', 'top') + elif isinstance(extrapolate, tuple): + extrapolate_options = [None, None] + for i in range(len(extrapolate)): + if isinstance(extrapolate[i], bool): + if extrapolate[i]: + extrapolate_options[i] = 'extrapolate' + else: + if i == 0: + extrapolate_options[i] = 'bottom' + else: + extrapolate_options[i] = 'top' + elif extrapolate[i] is None: + extrapolate_options[i] = np.nan + else: + extrapolate_options[i] = extrapolate[i] + extrapolate_options = tuple(extrapolate_options) + elif extrapolate is None: + extrapolate_options = ('bottom', 'top') + else: + extrapolate_options = (extrapolate, extrapolate) + + return extrapolate_options + + +def interpolate_vertical(self, new_levels, new_src_vertical=None, kind='linear', extrapolate_options=False, info=None, overwrite=False): """ Vertical interpolation. @@ -42,14 +105,35 @@ def interpolate_vertical(self, new_levels, new_src_vertical=None, kind='linear', the vertical data. kind : str Vertical methods type. - extrapolate : None, tuple, str - Extrapolate method (for non linear operations). + extrapolate_options : bool or tuple or None or number or NaN + If bool: + - If True, both extrapolation options are set to 'extrapolate'. + - If False, extrapolation options are set to ('bottom', 'top'). + If tuple: + - The first element represents the extrapolation option for the lower bound. + - The second element represents the extrapolation option for the upper bound. + - If any element is bool: + - If True, it represents 'extrapolate'. + - If False: + - If it's the first element, it represents 'bottom'. + - If it's the second element, it represents 'top'. + - If any element is None, it is replaced with numpy.nan. + - Other numeric values are kept as they are. + - If any element is NaN, it is kept as NaN. + If None: + - Both extrapolation options are set to (NaN, NaN). + If number: + - Both extrapolation options are set to the provided number. + If NaN: + - Both extrapolation options are set to NaN. info: None, bool Indicates if you want to print extra information. overwrite: bool Indicates if you want to compute the vertical interpolation in the same object or not. """ - + extrapolate_options = parse_extrapolate(extrapolate_options) + do_extrapolation = 'extrapolate' in extrapolate_options + if len(self.lev) == 1: raise RuntimeError("1D data cannot be vertically interpolated.") if not overwrite: @@ -126,7 +210,12 @@ def interpolate_vertical(self, new_levels, new_src_vertical=None, kind='linear', sys.stdout.flush() for j in range(ny): for i in range(nx): - curr_level_values = src_levels[t, :, j, i] + if len(src_levels.shape) == 1: + # To use 1D level information + curr_level_values = src_levels + else: + # To use 4D level data + curr_level_values = src_levels[t, :, j, i] try: # Check if all values are identical or masked if ((isinstance(curr_level_values, np.ndarray) and @@ -136,15 +225,24 @@ def interpolate_vertical(self, new_levels, new_src_vertical=None, kind='linear', kind = 'slinear' else: kind = kind # 'cubic' - if extrapolate is None: + + # Filtering filling values to extrapolation + fill_value = [np.nan, np.nan] + if 'bottom' in extrapolate_options: + if ascendant: + fill_value[0] = np.float64(self.variables[var_name]['data'][t, 0, j, i]) + else: + fill_value[0] = np.float64(self.variables[var_name]['data'][t, -1, j, i]) + else: + fill_value[0] = extrapolate_options[0] + if 'top' in extrapolate_options: if ascendant: - fill_value = (np.float64(self.variables[var_name]['data'][t, 0, j, i]), - np.float64(self.variables[var_name]['data'][t, -1, j, i])) + fill_value[1] = np.float64(self.variables[var_name]['data'][t, -1, j, i]) else: - fill_value = (np.float64(self.variables[var_name]['data'][t, -1, j, i]), - np.float64(self.variables[var_name]['data'][t, 0, j, i])) + fill_value[1] = np.float64(self.variables[var_name]['data'][t, 0, j, i]) else: - fill_value = extrapolate + fill_value[1] = extrapolate_options[1] + fill_value = tuple(fill_value) # We force the methods with float64 to avoid negative values # We don't know why the negatives appears with float34 @@ -155,24 +253,50 @@ def interpolate_vertical(self, new_levels, new_src_vertical=None, kind='linear', # 4D vertical component src_levels_aux = src_levels[t, :, j, i] - if kind == 'linear' and ascendant: + if kind == 'linear' and ascendant and not do_extrapolation: dst_data[t, :, j, i] = np.array( np.interp(new_levels, np.array(src_levels_aux, dtype=np.float64), - np.array(self.variables[var_name]['data'][t, :, j, i], dtype=np.float64)), + np.array(self.variables[var_name]['data'][t, :, j, i], dtype=np.float64), + left=fill_value[0], right=fill_value[1]), dtype=self.variables[var_name]['data'].dtype) else: - dst_data[t, :, j, i] = np.array( - interp1d(np.array(src_levels_aux, dtype=np.float64), - np.array(self.variables[var_name]['data'][t, :, j, i], dtype=np.float64), - kind=kind, - bounds_error=False, - fill_value=fill_value)(new_levels), - dtype=self.variables[var_name]['data'].dtype) + if not do_extrapolation: + dst_data[t, :, j, i] = np.array( + interp1d(np.array(src_levels_aux, dtype=np.float64), + np.array(self.variables[var_name]['data'][t, :, j, i], dtype=np.float64), + kind=kind, + bounds_error=False, + fill_value=fill_value)(new_levels), + dtype=self.variables[var_name]['data'].dtype) + else: + # If extrapolation first we need to extrapolate all (below & above) + dst_data[t, :, j, i] = np.array( + interp1d(np.array(src_levels_aux, dtype=np.float64), + np.array(self.variables[var_name]['data'][t, :, j, i], + dtype=np.float64), + kind=kind, + bounds_error=False, + fill_value='extrapolate')(new_levels), + dtype=self.variables[var_name]['data'].dtype) + # Check values below the lower vertical level + if fill_value[0] != 'extrapolate': + if ascendant: + idx_bellow = np.where(new_levels < src_levels_aux[0]) + else: + idx_bellow = np.where(new_levels > src_levels_aux[0]) + dst_data[t, idx_bellow, j, i] = fill_value[0] + # Check values above the upper vertical level + if fill_value[1] != 'extrapolate': + if ascendant: + idx_above = np.where(new_levels > src_levels_aux[-1]) + else: + idx_above = np.where(new_levels < src_levels_aux[-1]) + dst_data[t, idx_above, j, i] = fill_value[1] except Exception as e: print("time lat lon", t, j, i) print("***********************") - print("LEVELS", np.array(src_levels[t, :, j, i], dtype=np.float64)) + print("LEVELS", src_levels_aux) print("DATA", np.array(self.variables[var_name]['data'][t, :, j, i], dtype=np.float64)) print("METHOD", kind) print("FILL_VALUE", fill_value) @@ -191,13 +315,15 @@ def interpolate_vertical(self, new_levels, new_src_vertical=None, kind='linear', if flip: self.reverse_level_direction() new_lev_info['positive'] = self._lev['positive'] - for var_attr, attr_info in self.variables[self.vertical_var_name].items(): - if var_attr not in ['data', 'dimensions', 'crs', 'grid_mapping']: - new_lev_info[var_attr] = copy(attr_info) - self.set_levels(new_lev_info) - self.free_vars(self.vertical_var_name) - self.vertical_var_name = None + if self.vertical_var_name is not None: + for var_attr, attr_info in self.variables[self.vertical_var_name].items(): + if var_attr not in ['data', 'dimensions', 'crs', 'grid_mapping']: + new_lev_info[var_attr] = copy(attr_info) + self.free_vars(self.vertical_var_name) + self.vertical_var_name = None + + self.set_levels(new_lev_info) # Remove original file information self.__ini_path = None diff --git a/nes/nc_projections/default_nes.py b/nes/nc_projections/default_nes.py index 6c4ad49ee62109e4b79320cbdf8d2a1961f8176e..f81b8fe15abb4f0ec2baea3a2f6e5ea38e3ea605 100644 --- a/nes/nc_projections/default_nes.py +++ b/nes/nc_projections/default_nes.py @@ -3608,14 +3608,33 @@ class Nes(object): self : Nes Source Nes object. new_levels : List - New vertical levels. + List of new vertical levels. new_src_vertical : nes.Nes, str Nes object with the vertical information as variable or str with the path to the NetCDF file that contains the vertical data. kind : str Vertical methods type. - extrapolate : None, tuple, str - Extrapolate method (for non linear operations). + extrapolate : bool or tuple or None or number or NaN + If bool: + - If True, both extrapolation options are set to 'extrapolate'. + - If False, extrapolation options are set to ('bottom', 'top'). + If tuple: + - The first element represents the extrapolation option for the lower bound. + - The second element represents the extrapolation option for the upper bound. + - If any element is bool: + - If True, it represents 'extrapolate'. + - If False: + - If it's the first element, it represents 'bottom'. + - If it's the second element, it represents 'top'. + - If any element is None, it is replaced with numpy.nan. + - Other numeric values are kept as they are. + - If any element is NaN, it is kept as NaN. + If None: + - Both extrapolation options are set to (NaN, NaN). + If number: + - Both extrapolation options are set to the provided number. + If NaN: + - Both extrapolation options are set to NaN. info: None, bool Indicates if you want to print extra information. overwrite: bool @@ -3623,7 +3642,7 @@ class Nes(object): """ return vertical_interpolation.interpolate_vertical( - self, new_levels, new_src_vertical=new_src_vertical, kind=kind, extrapolate=extrapolate, info=info, + self, new_levels, new_src_vertical=new_src_vertical, kind=kind, extrapolate_options=extrapolate, info=info, overwrite=overwrite) def interpolate_horizontal(self, dst_grid, weight_matrix_path=None, kind='NearestNeighbour', n_neighbours=4, diff --git a/tests/3.1-test_vertical_interp.py b/tests/3.1-test_vertical_interp.py index a28567559832a7af81e38f61aa7785f17218dcfe..bbee5781a4b7fa7bf2d13f3347f5d08e46ab2e41 100644 --- a/tests/3.1-test_vertical_interp.py +++ b/tests/3.1-test_vertical_interp.py @@ -15,7 +15,7 @@ parallel_method = 'T' result_path = "Times_test_3.1_vertical_interp_{0}_{1:03d}.csv".format(parallel_method, size) result = pd.DataFrame(index=['read', 'calculate'], - columns=['3.1.1.Interp']) + columns=['3.1.1.Interp', '3.1.1.Exterp']) # ====================================================================================================================== # =============================================== VERTICAL INTERPOLATION ============================================= @@ -61,6 +61,50 @@ if rank == 0: print(result.loc[:, test_name]) sys.stdout.flush() +# ====================================================================================================================== +# =============================================== VERTICAL INTERPOLATION ============================================= +# ====================================================================================================================== + +test_name = '3.1.1.Exterp' +if rank == 0: + print(test_name) + +# READ +st_time = timeit.default_timer() + +# Original path: /esarchive/exp/monarch/a4dd/original_files/000/2022111512/MONARCH_d01_2022111512.nc +# Rotated grid from MONARCH +source_path = '/gpfs/projects/bsc32/models/NES_tutorial_data/MONARCH_d01_2022111512.nc' + +# Read source data +source_data = open_netcdf(path=source_path, info=True) + +# Select time and load variables +source_data.keep_vars(['O3', 'mid_layer_height_agl']) +source_data.load() + +comm.Barrier() +result.loc['read', test_name] = timeit.default_timer() - st_time + +# INTERPOLATE +st_time = timeit.default_timer() +source_data.vertical_var_name = 'mid_layer_height_agl' +level_list = [0., 50., 100., 250., 500., 750., 1000., 2000., 3000., 5000., 21000, 25000, 30000, 40000, 50000] +interp_nes = source_data.interpolate_vertical(level_list, info=True, kind='linear', extrapolate=True) +comm.Barrier() +result.loc['calculate', test_name] = timeit.default_timer() - st_time + +# WRITE +st_time = timeit.default_timer() +interp_nes.to_netcdf(test_name.replace(' ', '_') + "_{0:03d}.nc".format(size)) +comm.Barrier() +result.loc['write', test_name] = timeit.default_timer() - st_time + +comm.Barrier() +if rank == 0: + print(result.loc[:, test_name]) +sys.stdout.flush() + if rank == 0: result.to_csv(result_path) print("TEST PASSED SUCCESSFULLY!!!!!") diff --git a/tutorials/4.Interpolation/4.1.Vertical_Interpolation.ipynb b/tutorials/4.Interpolation/4.1.1.Vertical_Interpolation.ipynb similarity index 100% rename from tutorials/4.Interpolation/4.1.Vertical_Interpolation.ipynb rename to tutorials/4.Interpolation/4.1.1.Vertical_Interpolation.ipynb diff --git a/tutorials/4.Interpolation/4.1.2.Vertical_Interpolation_Extrapolation.ipynb b/tutorials/4.Interpolation/4.1.2.Vertical_Interpolation_Extrapolation.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..128ba78221724a008df2e10259a03b503558f20e --- /dev/null +++ b/tutorials/4.Interpolation/4.1.2.Vertical_Interpolation_Extrapolation.ipynb @@ -0,0 +1,302 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Vertical interpolation - Extrapolation options" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from nes import *\n", + "import matplotlib.pyplot as plt\n", + "from datetime import datetime\n", + "import numpy as np\n", + "\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "ascendant = True" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Creating base" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "var_name = 'var1'" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "lat_orig = 33.0\n", + "lon_orig = -15.0\n", + "inc_lat = 1.0\n", + "inc_lon = 1.0\n", + "n_lat = 13\n", + "n_lon = 25\n", + "\n", + "nessy = create_nes(info=False, projection='regular', \n", + " lat_orig=lat_orig, lon_orig=lon_orig, \n", + " inc_lat=inc_lat, inc_lon=inc_lon, \n", + " n_lat=n_lat, n_lon=n_lon)\n", + "\n", + "nessy.set_time([datetime(year=1986, month=11, day=29)])\n", + "nessy.set_levels({'data': (np.linspace(0, 2*np.pi, 10) * 10) + 100, 'units': 'm'})" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "if not ascendant:\n", + " nessy.lev['data'] = np.flip(nessy.lev['data'] )\n", + " print(nessy.lev)" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "nessy.variables[var_name] = {'data': None,\n", + " 'units': 'kg.m-2',\n", + " 'description': 'new variable from scratch',\n", + " 'short_name': 'newvar2',\n", + " }\n", + "\n", + "# Fill the data\n", + "shape = (len(nessy.time), len(nessy.lev['data']), nessy.lat['data'].shape[0], nessy.lon['data'].shape[-1])\n", + "nessy.variables[var_name]['data'] = np.zeros(shape)\n", + "\n", + "for y in range(nessy.lat['data'].shape[0]):\n", + " for x in range(nessy.lon['data'].shape[-1]):\n", + " nessy.variables[var_name]['data'][0,:,y,x] = np.sin(nessy.lev['data'] -100) +1" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deVyVZf7/8dcHERV3c02B475vsQhljcugZZlly5jkt22isqlmmmnGojJ/E5OzVzPaSGlqkVZTmNMqWlkWJGrlxlhuIGpKSS6ZAvL5/XEOZ1ARjsg5N8vn+XjwiHPd933OG+YePl73dd/XJaqKMcYYAxDkdABjjDE1hxUFY4wxXlYUjDHGeFlRMMYY42VFwRhjjJcVBWOMMV51qiiIyKUiskVEtorItHK2i4g87dm+XkQu8PVYP+dK8ORZLyKfisjgMtt2isgGEflCRNZUZy4fs40QkYOez/9CRB719dgAZHugTK6NInJCRNp4tvnt9yYi80Rkv4hsPMN2p86zynI5eZ5Vls3J86yybE6dZ2Ei8oGIZIvIJhG5r5x9qv9cU9U68QU0ALYB3YAQ4Eug3yn7jAPeAQSIBT7z9Vg/57oQaO35/rLSXJ7XO4G2Dv7ORgBvVuVYf2c7Zf/xwPsB+r1dAlwAbDzD9oCfZz7mcuQ88zGbI+eZL9kcPM86ARd4vm8OfBWIv2l1qacQA2xV1e2qWggsBiacss8EYKG6ZQKtRKSTj8f6LZeqfqqqBZ6XmUCXavrsc87mp2P98f43AIuq8fPPSFU/Ag5UsIsT51mluRw8z3z5nZ2Jv8+zs80WyPNsr6qu83x/GMgGOp+yW7Wfa3WpKHQGdpV5ncfpv8Az7ePLsf7MVdZtuCt/KQWWichaEUmspkxnmy1ORL4UkXdEpP9ZHuvvbIhIKHAp8FqZZn/+3irjxHl2tgJ5nvnKifPMZ06eZyLiAoYCn52yqdrPteCqhqyBpJy2U+fwONM+vhxbVT6/t4iMxP1/1uFlmi9S1T0i0h5IF5H/ev5lE6hs64AIVT0iIuOAJUBPH4/1d7ZS44FPVLXsv/b8+XurjBPnmc8cOM984dR5djYcOc9EpBnuQvRLVT106uZyDjmnc60u9RTygLAyr7sAe3zcx5dj/ZkLERkEPAdMUNXvSttVdY/nv/uBNNzdwupSaTZVPaSqRzzfvw00FJG2vhzr72xlTOKULr2ff2+VceI884lD51mlHDzPzkbAzzMRaYi7IKSq6uvl7FL955o/Bkic+MLd69kOdOV/Ayv9T9nnck4elFnt67F+zhUObAUuPKW9KdC8zPefApcG+HfWERDP9zFAruf357ff2dn8bwK0xH09uGmgfm+e93Vx5kHTgJ9nPuZy5DzzMZsj55kv2Zw6zzw//0LgyQr2qfZzrc5cPlLVYhH5BfAe7pH3eaq6SUTu9Gz/F/A27tH6rcBR4JaKjg1grkeB84DZIgJQrKpRQAcgzdMWDLykqu9WR66zyHYtcJeIFAM/ApPUfdb57Xd2FtkArgaWqeoPZQ736+9NRBbhvlumrYjkAdOBhmVyBfw88zGXI+eZj9kcOc98zAYOnGfARcAUYIOIfOFpewh3cffbuVZamY0xxpg6NaZgjDHmHFlRMMYY42VFwRhjjJcVBWOMMV71rig4+LRmpWpqtpqaCyxbVdTUXGDZqqK6c9W7ogDUyP9hPWpqtpqaCyxbVdTUXGDZqsKKgjHGGP+o1c8ptG3bVl0u11kdk5+fT7t27fwT6BzV1Gw1NRdYtqqoqbnAslVFVXKtXbv2W1Ut96Ba/USzy+VizZpqXw/EGGPqNBHJOdM2u3xkjDHGy4qCMcYYLysKxhhjvKwoGGOM8bKiYIwxxsuKgjHVLDU1FZfLRVBQEC6Xi9TUVKcjGeOzWn1LqjE1TWpqKomJiRw9ehSAnJwcEhPdD5wmJCQ4Gc0Yn9Tqh9eioqLUnlMwNUl4eDi7du06rT0iIoKdO3cGPpAx5RCRtZ5V905jPQVjqkhV+eqrr8jIyCAzM5OMjIxyCwJAbm5ugNMZUzVWFIzx0aFDh1i9ejXnnXceQ4cO5auvvqJPnz4AtGzZkmHDhpGTk8PBgwdPOzY8PDzQcY2pEhtoNuYMVJV58+Zx++23M3DgQFq1akV8fDyzZs0CoGfPnsybN4+NGzdy4MAB3nvvPWbNmkVoaOhJ7xMaGkpycrITP4IxZ83GFIwBvv/+ez777DMyMzMRER599FEAevXqRX5+PrGxscTFxREbG0tMTAytWrU643ulpqaSlJREbm4u4eHhJCcn2yCzqVEqGlOwomDqHVVFRACYOXMmCxcuJDs7GwARYeTIkaxYsQKAffv20a5dO4KCrFNt6g4baDb12oEDB/jss8+8A8Lr168nNzeXkJAQjh8/Trdu3Zg8eTJxcXFER0fTokUL77EdOnRwMLkxgee3oiAi84ArgP2qOqBM+z3AL4Bi4C1V/a2n/UHgNuAEcK+qvuevbKbuOnHiBJs2baJbt240a9aM2bNnc/fddwMQFBTEoEGDmDhxIkeOHKFNmzZMnz7d4cTG1Cz+7CnMB/4JLCxtEJGRwARgkKoeF5H2nvZ+wCSgP3A+sFxEeqnqCT/mM3XAkSNH+PDDD723hK5evZojR47w1ltvMW7cOC666CL+8Ic/EBcXR1RUFM2aNXM6sjE1mt+Kgqp+JCKuU5rvAmaq6nHPPvs97ROAxZ72HSKyFYgBMvyVz9Q+xcXFbNiwgYyMDAYNGsTw4cPZsWMH48ePp0GDBgwePJibbrqJ2NhYoqLcl0sHDx7M4MGDHU5uTO0R6DGFXsDFIpIMHAN+o6pZQGcgs8x+eZ6204hIIp6Fqu3e77qvuLiYRx55hIyMDLKysrzTR9x///0MHz6cfv368eGHHxIVFUXTpk0dTmtM7RfoohAMtAZigWjgFRHpBkg5+5Z7W5SqpgAp4L77yE85TYAVFRWxfv1672Bw27ZtefLJJwkODua1116jVatW/PznP/feGhoREQFAgwYN+MlPfuJwemPqjkAXhTzgdXXfB7taREqAtp72sDL7dQH2BDibCaCCggJat24NwB133MHChQs5duwYAJ06deKqq67y7pudnU2DBg0cyWlMfRPoorAEGAV8KCK9gBDgW2Ap8JKI/A33QHNPYHWAsxk/KSws5IsvvjhpjqD8/HwOHjxIcHAwffv25c477/Q+HBYWFuZ9jgCwgmBMAPnzltRFwAigrYjkAdOBecA8EdkIFAI3eXoNm0TkFWAz7ltV77Y7j2qv3bt3k5mZyZgxY2jevDlPPPEEjz32GABdunQhLi6OuLg4CgsLCQ4O5pe//KWzgY0xXvZEszlne/bs4eWXX/b2BEpnCl22bBnx8fFkZ2ezadMmYmNj6dKli8NpjTH2RLOpFqpKXl4eGRkZZGRkMG7cOOLj4/nmm2+4//77iYiI4KKLLvJeBhoyZAgAffv2pW/fvg6nN8b4woqCOaPSOYJ+/PFHpkyZQkZGBnv2uMf/GzduTHh4OPHx8QwaNIg9e/bQqVMnhxMbY86VFQUDuAtATk6OtxeQmZlJv379mD9/Po0bN2bv3r2MGDHCe0vo4MGDadiwIQDBwcFWEIypI6wo1FNHjx5l27ZtDBw4EID4+HjvzKChoaFER0czaNAgwD1z6CeffOJYVmNM4Nh8wHVAamoqLpeLoKAgXC4Xqampp+2Tl5dHamoq99xzD1FRUbRs2ZLhw4dTUlICwI033sisWbNYt24dBw8e5MMPP+T+++8P9I9ijHGY9RRqudTUVBITE73TP+Tk5HD77beTnZ1N06ZNueeee2jWrBn/+te/SE5OplmzZsTExPDb3/6WuLg4SkpKCAoK4uabb3b2BzHG1Ah2S2ot53K5yMnJOeP2Tz/9lLi4OO/awf3797eHwYyp5+yW1DrqyJEj5ObmlrtNRMjPz+e8884D8M4VZIwxFbExhVpo8+bN3HPPPZx//vmnLRJfKjw83FsQjDHGV1YUapE333yTUaNG0b9/f1JSUpgwYQIPPPDAaYUhNDSU5ORkh1IaY2ozu3xUw+3du5eOHTsiIqxYsYLt27fzxBNPcNttt9GuXTsAevToQVJSErm5uYSHh5OcnExCQoLDyY0xtZENNNdAqsrKlSuZNWsWaWlpLFu2jFGjRnH48GFCQ0NtoNgYc05soLmWOHbsGHPnzmX27Nls3ryZ1q1b88tf/pLu3bsD0Lx5c4cTGmPqOisKNUDpgjNBQUE8/vjjdO7cmXnz5jFp0iSaNGnidDxjTD1iRcEhhYWFpKWlMWvWLHbu3Mn27dsJCQnh888/p2PHjk7HM8bUU1YUAmzPnj0888wzPPvss+zbt4+uXbtyzz33UFxcTHBwsBUEY4yjrCgEgKpy7NgxmjRpwpdffklycjKXX345U6dOZezYsQQF2Z3Bxpiawf4a+dH333/P008/Td++fZkxYwYAY8eOZfv27fznP//hsssus4JgjKlR7C+SH3z55ZckJibSuXNn7rvvPlq1akV0dDSAdyZTY4ypiezyUTUpKiryLjrzpz/9ibS0NCZPnsxdd91FZGSkw+mMMcY31lM4Rzk5OTz00EN07tyZDRs2ADBz5kx2797Nc889ZwXBGFOrWE+hCkpKSkhPT2f27Nm8+eabAIwfPx4RASAsLMzJeMYYU2V+6ymIyDwR2S8iG8u0PSYiu0XkC8/XuDLbHhSRrSKyRUTG+ivXuShdpezw4cNcc801ZGRkMG3aNLZv386SJUsYMGCAwwmNMebc+LOnMB/4J7DwlPa/q+pfyjaISD9gEtAfOB9YLiK9VPWEH/P5bO3atcyePZvs7Gw++eQTWrZsyfvvv8/gwYNp1KiR0/GMMaba+K2noKofAQd83H0CsFhVj6vqDmArEOOvbL44duwYCxcuZNiwYURFRbF48WIGDhzIsWPHAIiJibGCYIypc5wYaP6FiKz3XF5q7WnrDOwqs0+epy3gSmeNffXVV7nppps4dOgQTz/9NHv27GHOnDk2F5Expk4LdFF4BugODAH2An/1tEs5+5Y7p7eIJIrIGhFZk5+fX6UQqampuFwu7zMDL7zwAm+//TZXXHEFTz31FADXXnsty5cv965y1rJlyyp9ljHG1CYBvftIVfeVfi8izwJvel7mAWVv2ekC7DnDe6QAKeBeT+FsM6SmppKYmMjRo0cB9y2lN910E6pKhw4dGDfOPfbdpEkTRo8efbZvb4wxtVpAi4KIdFLVvZ6XVwOldyYtBV4Skb/hHmjuCaz2R4akpCRvQSilqrRt25bc3FxCQkL88bHGGFMr+K0oiMgiYATQVkTygOnACBEZgvvS0E7gDgBV3SQirwCbgWLgbn/deZSbm1tu+3fffWcFwRhT7/mtKKjqDeU0z61g/2TA76vNh4eHk5OTU267McbUd/Vumovk5GRCQ0NPamvYsCHJyX6vR8YYU+PVu6KQkJBASkoKERERiAihoaEUFRXRtGlTp6MZY4zj6l1RAHdh2LlzJyUlJXz77bfExMRw44038sUXXzgdzRhjHFUvi0JZTZo0YcmSJfTr14+ioiKn4xhjjKNsllSgU6dOfPbZZ95ZTktKSmxFNGNMvWR/+TxKC0JSUhL/93//553uwhhj6hMrCqdo1qwZqampPPHEE05HMcaYgLPLR6eYNm0amzdvJikpiT59+jBx4kSnIxljTMBYT+EUIsKzzz5LbGwsU6ZMYd26dU5HMsaYgKm0KIhInIjM8kx3nS8iuSLytojcLSJ1curQxo0bs2TJEjp37sy2bducjmOMMQFT4eUjEXkH92ylb+CegmI/0BjoBYwE3hCRv6nqUn8HDbQOHTqwceNGmw/JGFOvVDamMEVVvz2l7QiwzvP1VxFp65dkNUBpQXjppZdYtmwZzz//vPcuJWOMqYsqvHx0akEQkRYi0qb0q7x96qK8vDwWLFjA73//e6ejGGOMX/l095GI3AH8P+BH/rcimgLd/JSrRnnggQfYvHkz06dPp0+fPlx//fVORzLGGL/w9ZbU3wD960OvoDwiwpw5c9i6dSs33XQTXbt2JTo62ulYxhhT7Xy9JXUbcLTSveqwRo0akZaWRseOHUlPT3c6jjHG+IWvPYUHgU9F5DPgeGmjqt7rl1Q1VLt27fj8889p1aqV01GMMcYvfO0pzAHeBzKBtWW+6p3SgrBmzRri4+OJiIggKCgIl8tFamqqw+mMMebc+NpTKFbV+/2apJZ58sknWb58ufd1Tk4OiYmJgHu9BmOMqY187Sl8ICKJItLp1FtS66uPP/74tLajR4+SlJTkQBpjjKkevvYUJnv++2CZtnpzS2p5du3aVW57bm4uW7ZsoVGjRrhcrsCGMsaYc+RTT0FVu5bzVW8LAkB4ePgZ26dPn07Xrl3p1asXd999N2+88QaHDh0KcEJjjDl7FRYFERleyfYWIjKgeiPVDsnJyYSGhp7UFhoaSnJyMjNmzOCpp56id+/eLFiwgKuuuopLLrnEu192djbFxcWBjmyMMZWq7PLRNSLyJ+Bd3Hcb5eOeEK8H7gnxIoBfl3egiMwDrgD2q+qAU7b9Bvgz0K70gTgReRC4DTgB3Kuq71X1hwqE0sHkpKQkcnNzCQ8PJzk52dveu3dv7r33XgoLC8nMzOSHH34AoKioiGHDhiEijBw5kjFjxhAfH0+PHj1sXiVjjOOksmUnRaQ1cC1wEdAJ91QX2cBbqrqqguMuwT153sKyRUFEwoDngD5ApKp+KyL9gEVADHA+sBzopaonKsoWFRWla9asqfSHrEkKCwtZunQp6enpLFu2jJ07dwLwxBNPMG3aNAoLCzly5Aht2tTrcXxjjB+JyFpVjSpvW6UDzapaADzr+fKZqn4kIq5yNv0d+C3u6bhLTQAWq+pxYIeIbMVdIDLO5jNrg5CQEK699lquvfZaVJVt27axbNkyhg93X6lbuXIlY8eOJTo6mvj4eOLj44mLi7MpvI0xARHQlddE5Epgt6p+ecqmzkDZ23nyPG3lvUeiiKwRkTX5+fl+ShoYIkKPHj2YOnUqgwYNAqB79+48+uijBAcHM3PmTEaMGEGbNm3YsmULAD/++COV9e6MMaaqAlYURCQUSAIeLW9zOW3l/uVT1RRVjVLVqHbt2lVnxBqhW7duPPbYY3zyySd89913pKWl8fOf/5wePXoA8Lvf/Y7w8HBuvfVWFi9eTG0vjMaYmsXX5xSqQ3egK/ClZ0C1C7BORGJw9wzCyuzbBfeKb/Vay5Ytueqqq7jqqqu8bZdccgl79uwhLS2N559/HoDLLruMt99+G4CSkhKCgmzpbWNM1fi6nsJ1wLuqelhEHgYuAB5XVZ9XtVfVDUD7Mu+5E4jyDDQvBV4Skb/hHmjuCaz2/ceoP0rHI06cOMGaNWtIT0+nQYMGAKgq/fv3x+VyER8fz5gxY+jfv7/d1WSM8ZmvPYVHVPVVz3MLY4G/AM8Aw850gIgsAkYAbUUkD5iuqnPL21dVN4nIK8BmoBi4u7I7j+q7Bg0aMGzYMIYN+9//BMePH2fMmDGkp6fz61+77xTu1KkTycnJ3HLLLU5FNcbUIr4WhdI/0JcDz6jqGyLyWEUHqOoNlWx3nfI6GUj2MY8pR+PGjXnqqacA9zQcy5cvZ9myZbRv7+6grV+/nilTpnh7ERdffDFNmjRxMrIxpoap9DkFABF5E9gN/BSIxP2swmpVHezfeBWrjc8pOGn16tU89NBDfPzxxxQWFtKoUSMuvvhi5syZQ7du9XrWEmPqlYqeU/B1RPJ64D3gUlX9HmgDPFBN+UyAxMTEsHz5cgoKCnj33Xe5++67+fbbbym9i+vpp59m8uTJPP/88+Tl5Tmc1hjjhAovH50yPfaHZdqOA/ZP9FoqNDSUsWPHMnbs2JPajxw5wvvvv8+iRYsA6Nu3LxMmTOCJJ55wIqYxxgEVXj4SkR24nxco9zkCp2dKtctH1U9V2bBhg3cajtDQUNLS0gC48847CQsLY8yYMVxwwQXeu56MMbVLRZePfBpTqKmsKPifqiIiFBYWEhcXx7p17ruQW7duzejRo0lMTCQ+Pt7hlMaYs3HOYwridqOIPOJ5He556MzUcaXPOISEhLB27Vr27dtHamoqEyZM4NNPP+Xrr78GYO/evbZ2hDF1gK93Hz0DlACjVLWvZ+bUZaoa7e+AFbGegrNUleLiYho2bMjy5cuZMGECR48epUGDBsTGxhIfH88dd9xBx44dnY5qjCmjOu4+GqaqdwPHwDtzqk3bWc+JCA0bNgTgpz/9KQcOHOCDDz7gd7/7HcePH2fGjBkUFRUBsGzZMp555hm2bdt20nukpqbicrkICgrC5XKRmpoa8J/DGPM/vvYUPgMuBLJU9QIRaYe7pzDU3wErYj2Fmq2goIDWrVsDcMcdd5CSkgJA165dGTNmDE2aNCElJYWjR496jwkNDSUlJcW7WJExpvqd80CziCQAP8M959EC3IvuPKyqr1Zn0LNlRaH2UFW2bt3qvavp/fff59ixY96eRFnh4eHk5OQ4kNKY+qFa7j4SkT7AaNy3p65Q1ezqi1g1VhRqr6KiIho1anTGtSEuu+wybrzxRiZPnhzgZMbUfee08prnDZ4CXlbVWdWazNRbDRs2PGOPoGnTpuzatYutW7cCcOjQIQYPHkxkZCTR0dFER0cTFRVFixYtAh3bmDrP14HmdcDDIrJVRP4sIuVWGGPORnJyMqGhoSe1hYaGMmfOHDZs2MAjjzwCuIvCsGHD+Pzzz5k2bRqjR4+mVatWvPDCCwAcOHCAzMxMjh07FvCfwZi6xqeegqouABZ4pri4BvijiISrak+/pjN1WulgclJSErm5uYSHh5OcnOxtL31GokuXLixevBiAb7/9ljVr1pCVlUV0tPuO6GXLlnHDDTcQHBzMoEGDvL2JiRMnege6jTG+Oasnmj0PrP0MuArYrKrj/RXMFzamYMBdKD7++GNWr15NVlYWa9as4eDBg+zcuZOIiAhef/11Vq1a5S0W3bt3t4WHTL1WHXcf/RGYCGwDXgbSPLOlOsqKgilPSUkJW7dupWfPnogIM2bMYObMmd7LS23atGHYsGG8+eabBAUF8eOPP9q6EqZeqY6icCfwb1X9trrDnQsrCsZXRUVFbNq0iaysLLKysigoKODVV913VI8ZM4ZNmzYRExNz0kC2XXoydVV1FIUgYDLQTVX/n4iEAx1V1dF1lK0omOrw3HPP8cEHH5CVleWdy2nUqFGsWLECgIULF9K9e3eGDh162sC4MbVRdRQFm/vI1AsFBQWsXbuW4OBgRowYwY8//kjz5s05ceIEDRo0oH///kRHRzN58mRGjRrldFxjquScn1PAPffRBSLyObjnPhIRm/vI1DmtW7fmpz/9qfd1kyZN2LVrl/eyU1ZWFq+//jr9+/dn1KhR7N69m2uvvZbo6Gjv5aeePXsSFOTr3d7G1Cy+FoUiEWmAe8EdPHMflfgtlTE1SKdOnbjyyiu58sorAfeUHaXTcxQUFNCwYUPmzp3LP/7xDwBatmzJv//9b+8kgT/88ANdunSxO55MreBrUXgaSAPai0gy7rmPHvFbKmNqMBEhJMTdUR4wYAAfffQRxcXFZGdne3sTPXr0AODll19m6tSpdOjQwTuIHR0dzciRI2ncuLGTP4Yx5fLb3EciMg+4AtivqgM8bb8HJuDuZewHblbVPZ5tDwK3ASeAe1X1vcoy2ZiCqem2bt3Ku+++6y0W//3vf1FVCgoKaNWqFa+99ho7duwgJiaGCy64gGbNmjkd2dQDflmOU0RyVTW8gu2XAEeAhWWKQgtVPeT5/l6gn6reKSL9gEVADHA+sBzopaonKspgRcHUNocOHWLjxo1ceOGFANxyyy3Mnz8fgKCgIPr27cuIESP45z//CfxvOVRjqlN1DDSX+74VbVTVj0TEdUpb2XUam+IZo8Dde1isqseBHSKyFXeByDiHfMbUOC1atPAWBIDnn3+eP/7xj6xZs8b7RPbu3bu92y+55BKOHz9+0kB27969adCggRPxTT1wLkWhSl0Mz5jE/wEHgZGe5s5AZpnd8jxtxtR57du3Z9y4cYwbN+60bSNGjGDVqlW88MILzJ49G4BJkyaxaNEiAJYuXcqgQYOIiIiwHoWpFhUWBRG5/0ybgCpd/FTVJCDJM4bwC2A65fc6yi06IpIIJIJ7MRZj6rLf//73gHvqji1btpCVlUWnTp0A2LdvHxMmTACgXbt23kHsa665hoEDBzqW2dRulfUUmlew7alz/OyXgLdwF4U8IKzMti7AnvIOUtUUIAXcYwrnmMGYWqF0vKFv377etvPOO++ky05ZWVm88847uFwuBg4cyJYtW3j44Ye9xSIyMtLWoDCVqrAoqOqM6vwwEempql97Xl4J/Nfz/VLgJRH5G+6B5p6Ao1NoGFPTBQcHExkZSWRkJHfddRcAR44c8V5G2rNnD2vXruXf//434L6Vtnfv3ixevJjBgwdz+PBhGjZsaLfGmpP47bFLEVmEe6C4t4jkichtwEwR2Sgi64ExwH0AqroJeAXYDLwL3F3ZnUfGmNM1a9aMpk2bAjBy5Ei2b99Ofn4+77zzDjNmzKBnz56cf/75AMyaNYsWLVoQGRnJnXfeybx589iwYQMlJac/l5qamorL5SIoKAiXy0VqampAfy4TOFW+JbUmsFtSjam6zMxM3njjjZPWoAgJCeHw4cOEhITwxhtvcPToUb755hsefvhhjh496j02NDSUlJQU74JIpnbxy3MKNYEVBWOqR0lJCV9//TVbt27l8ssvB9w9jQ8//PCMx0RERLBz587ABDTVqspFoYK7jwBQ1b+dY7ZzYkXBGP8pKipi48aNXHDBBWfcZ+7cucTHxxMWFnbGfUzNU1FRqGxMoXklX8aYOqphw4YMHTqUiIiIcrcHBQVx2223ER4eTt++fbn33ntZu3ZtgFOa6hbQu4+MMbVPcnIyiYmJ5Y4pDBo0iPT0dNLT03nuuecYMmQIkZGR7Nq1iwULFjBmzBgiIyPtCexaxNdFdhrjnqyuP+C9f01Vb/VftMrZ5SNjAiM1NZWkpCRyc3MJDw8nOTn5tEHmY8eOocpta+8AABZaSURBVKo0adKEV155hZ/97GeAe42K0aNHEx8fz6RJk+xZiRqgOlZeexX3MwWTgf8HJADZqnpfdQY9W1YUjKm58vPzWbFiBenp6Sxbtoy8vDz27dtH+/btef/99zl48CAjR46kVatWTketd6qjKHyuqkNFZL2qDhKRhsB7quroeoRWFIypHVSVHTt20K1bNwCuvvpqlixZQoMGDYiJiSE+Pp5LL72UuLg4h5PWD+cy0FyqyPPf70VkANAScFVDNmNMPSAi3oIA7sWHVq5cybRp0zhx4gSPP/44v/nNb7zb09LS+Prrr6nNt8zXVr72FH4OvAYMBObjngzvUVX9l1/TVcJ6CsbUDQcOHGDfvn307duXH374gdatW1NUVERERARjxowhPj6e0aNH06ZNG6ej1gnn3FNQ1edUtUBVP1LVbqra3umCYIypO9q0aeOd7C80NJTNmzcza9Yshg4dyssvv8z1119PSkoKAIcPH2blypUUFhY6GbnO8mk9BRH5A/AnVf3e87o18GtVfdif4Ywx9Y+I0KNHD3r06MHUqVMpLi5m9erV3qny09PTueaaa2jatCkjRozw9iT69Olja0pUg7MaaD6lbZ2qnvlRxwCwy0fG1D+HDh3i/fff997VtHXrVgCys7Pp06cPOTk5hIaG0q5dO4eT1lzVsRxnAxFp5FkuExFpAjSqroDGGOOrFi1acNVVV3HVVVcBsGPHDlauXEnv3r0BmD59OgsWLGDo0KHeXsRFF11kU4T7yNeewm9xr3/wPO4V0W4Flqrqn/wbr2LWUzDGnOqLL77grbfeIj09nU8//ZSioiIGDhzI+vXrAcjLy6Nz5871+lJTtcySKiKXAaNxL525TFXfq76IVWNFwRhTkSNHjrBy5UqOHj3Kddddx4kTJ2jfvj0hISHEx8d7vzp27Oh01ICyqbONMQYoLCwkNTXVO1/Tt99+C8DMmTP53e9+x4kTJzh+/DihoaEOJ/Wvc5k6e5WqDheRw7gvG3k3Aaqqjk5iYkXBGFNVJSUlfPnllyxbtoyRI0cSExPDqlWrGD16NMOHD/eORwwZMoSgIL8tUukI6ykYY4wPvvrqK+bMmUN6ejobNmwAoG3btnz00Uf07duXEydO1IkZX8/54TURecGXNmOMqc169erFX//6V9avX8+ePXtYuHAh48ePp3v37gBMmzaNfv36cd999/HWW29x5MgRhxNXP1/7RP3LvhCRYCCy+uMYY0zN0KlTJ6ZMmcK8efMICQkBYMiQIYSHh5OSksIVV1xBmzZtuPrqqx1OWr0qLAoi8qBnPGGQiBzyfB0G9gFvBCShMcbUEAkJCbz77rsUFBSwfPlyfvWrX3mfjwC48MILue6663j22WdPW786NTUVl8tFUFAQLpeL1NTUAKf3TaVjCiISBDzn9II65bExBWNMTVFYWMgdd9xBeno6u3fvBqBnz548/PDDNGjQ4Iyr1526WFEgnNOYgqqWAIOrPZUxxtQhISEhPP/88+zatYvNmzfz1FNP0atXLxo3bkxSUtJJBQHg6NGjJCUlOZT2zHx9onkWMF9Vs3x+Y5F5wBXAflUd4Gn7MzAeKAS2AbeUmWTvQdxLfp4A7vXl4TjrKRhjaoOgoKBy14YQEUpKSgKepzoW2RkJZIrINhFZLyIbRGR9JcfMBy49pS0dGKCqg4CvgAc9AfsBk3APaF8KzBaR2n/flzHGgHeGV1/bneRrUbgM6AaMwv0v/Ss8/z0jVf0IOHBK2zJVLfa8zAS6eL6fACxW1eOqugPYCsT4mM0YY2q05OTk056SDg4OJjk52aFEZ+brIjs5QBgwyvP9UV+PrcCtwDue7zsDu8psy/O0nUZEEkVkjYisyc/PP8cIxhjjfwkJCaSkpBAREYGI0KJFC/r06cO1117rdLTT+Prw2nTgd3gu9wANgRer+qEikgQUA6X3ZJU3XWG5gx2qmqKqUaoaZfOlG2Nqi4SEBHbu3ElJSQkFBQV8/vnnNGpU81Yg8PVf+1fjnjr7BwBV3QM0r8oHishNuC8/Jej/Rl7ycPdESnUB9lTl/Y0xpqYLCgoiODiYb775hrFjx7J582anI3n5WhQKPX/AFUBEmlblw0TkUtw9jitVtez9WUuBSSLSSES6Aj2B1VX5DGOMqS0KCwv58ssvGT9+vHfGVqf5WhReEZE5QCsRuR1YDjxb0QEisgjIAHqLSJ6I3Ab8E3cPI11EvhCRfwGo6ibgFWAz8C5wt6qeqNJPZIwxtUR4eDhLlixh9+7dXHPNNRQWFjod6awW2YkHxuC+/v+eqqb7M5gv7DkFY0xd8NJLL5GQkMBtt93Gs88+6/dV4aq8RrOI/BN4SVU/9RQBxwuBMcbUNZMnTyY7O5uXX36ZgoIC2rRp41iWyi4ffQ38VUR2isgfRWRIIEIZY0x9M2PGDLKyshwtCFBJUVDVp1Q1DvgJ7gfRnheRbBF5VER6BSShMcbUA0FBQbRs2ZLjx49z1113sXHjRmdy+LKTquao6h9VdSgwGfctqtl+TWaMMfXQgQMHeOONNxg/fjz79+8P+Of7+vBaQxEZLyKpuJ9C/gq4xq/JjDGmHurUqRNLly7lm2++YeLEiRw/fjygn1/ZIjvxntlO84BE4G2gu6r+TFWXBCKgMcbUN1FRUSxYsIBPPvmEO++8s9wZVv2lsp7CQ7ifNeirquNVNVVVfwhALmOMqdeuv/56ZsyYQVpa2mmruPmTz88p1ET2nIIxpi5TVfLy8ggLC6t857NQHespGGOMCTARISwsDFXlySefZP36ypaxOXdWFIwxpob7/vvv+ctf/sL48ePZt2+fXz/LioIxxtRwrVu3ZunSpeTn5zN8+HDCw8MJCgrC5XKRmppa+RucBSsKxhhTC1xwwQXcfvvtbN26lV27dqGq5OTkkJiYWK2FwYqCMcbUEm+88cZpbUePHiUpKanaPqPCCfGMMcY47/Dhw7z44ovk5OSUuz03N7faPsuKgjHG1FCbNm1i9uzZLFy4kCNHjtCoUaNyn3AODw+vts+0y0fGGFMDzZw5kwEDBjB37lwmTpxIZmYmc+fOJTQ09KT9QkNDSU5OrrbPtZ6CMcbUALt37yYlJYWJEycyePBgLr30UoKCgrj11ltp27YtAMOGDQMgKSmJ3NxcwsPDSU5OJiEhodpyWFEwxhiHqCoffPABs2fPZsmSJZSUlNC6dWsGDx7MkCFDGDLk9CVsEhISqrUInMqKgjHGOEBVufDCC8nMzKRNmzb86le/4s4776R79+6O5rIxBWOMCZD169fz6KOPoqqICJMmTWL+/Pnk5eXx5z//2fGCADYhnjHG+FVhYSGvvfYas2fPZtWqVTRu3JgvvviC3r17O5bJJsQzxhgHbNiwgbCwMCZPnszevXv5y1/+wu7dux0tCJXxW1EQkXkisl9ENpZpu05ENolIiYhEnbL/gyKyVUS2iMhYf+Uyxhh/KSkpIT09nddffx2A3r17Ex8fzzvvvMNXX33Fr3/9a9q0aeNwyor5c6B5PvBPYGGZto3ARGBO2R1FpB8wCegPnA8sF5FeqnrCj/mMMaZaFBQUsGDBAp555hm++uorIiMjmThxIiEhIbz44otOxzsrfuspqOpHwIFT2rJVdUs5u08AFqvqcVXdAWwFYvyVzRhjqsusWbPo3Lkzv/rVrzjvvPN44YUXWLVqldOxqqymjCl0BnaVeZ3naTuNiCSKyBoRWZOfnx+QcMYYU+rYsWO8+OKL3vmGevbsSUJCAuvWrePTTz/lxhtvpHHjxg6nrLqaUhSknLZyb4tS1RRVjVLVqHbt2vk5ljHGuO3cuZMHH3yQsLAwpkyZwksvvQTAmDFjePbZZxk6dKjDCatHTXl4LQ8ouwhpF2CPQ1mMMcarpKSE6667jrS0NESEK6+8kqlTpzJ69Gino/lFTekpLAUmiUgjEekK9ARWO5zJGFNPHThwgEWLFgEQFBRE586deeihh9i5cydpaWnEx8cTFFRT/nxWL7/1FERkETACaCsiecB03APP/wDaAW+JyBeqOlZVN4nIK8BmoBi42+48MsYE2po1a5g1axaLFy/m2LFjDBs2jG7duvH00087HS1g/FYUVPWGM2xKO8P+yUD1zf9qjDE+2rx5MzfffDNZWVk0bdqUm2++mbvuuotu3bo5HS3gasqYgjHGBNT27dvJz89n2LBhnH/++agq//jHP5gyZQotW7Z0Op5j6uZFMWOMAVJTU3G5XAQFBeFyuXjhhRd46623GDduHD169OCee+4BoFWrVmRlZfGLX/yiXhcEsJ6CMaaOSk1NJTExkaNHjwKQk5PDTTfdhKrSqVMnHn30UW6//XaHU9Y8VhSMMXXOkSNHuP/++70FoZSq0rZtW3JycmjYsKFD6Wo2u3xkjKkT8vLyuOuuuxgyZAgtW7Zk//795e733XffWUGogBUFY0ytcujQIVasWMHjjz/O5Zdfzr/+9S8AGjZsyEsvvUT79u15+OGHad++fbnHh4eHBzJurWOXj4wxNVZJSQkFBQWcd955lJSUEBMTw7p16yhdHKxfv36UlJQA0KFDBwoKCrwPlfXq1eukMQWA0NBQkpPtzveKWFEwxtQY33//PatXryYjI4OMjAw+++wzBg0axMqVKwkKCuKiiy5i/PjxxMbGMmzYMFq1anXS8WWfMi5d3D4pKYnc3FzCw8NJTk7266L3dYEtx2mMcURJSQnZ2dls2rSJ66+/HoDLL7+ct99+GxFhwIABxMbGMmLECCZPnuxw2rqlouU4rSgYYwLmyy+/5PXXX/f2Ag4dOgS4F6lp1aoVH3/8McePHycmJoYWLVo4nLbuqqgo2OUjY0y1O3HiBJs2bSIzM5OMjAymT5+Oy+UiIyODxx9/nIEDBzJ58mRiY2OJi4vzPjB28cUXO5zcWE/BGHPOVBURYdOmTdx3332sXr2aw4cPA9C2bVteeeUVRo4cyeHDhxERmjVr5nDi+s16CsaYalNcXMzGjRvJyMggMzOTzMxMpk6dyn333UfLli05cOAAU6ZMIS4ujtjYWLp3746Iex2t5s2bO5zeVMaKgjGmQvn5+Rw4cIDevXtTWFhI+/btOXjwIADt27cnLi6Orl27AtClSxfWrVvnZFxzjqwoGGNOsn79ej7++GNvT2Dbtm2MGjWKFStWEBISwrRp04iIiCA2NhaXy+XtBZi6wcYUjKnHvvnmGzIyMtixYwf3338/APHx8SxfvpyOHTsSFxdHXFwcF198MbGxsQ6nNdXFbkk1xni99957zJ8/n8zMTHbu3AlAkyZN+O6772jSpAkbN26kefPmhIeHWy+gjrKBZmPqoT179ngvAWVkZLBo0SLCwsLYsmULq1atIjY2lnvuuYe4uDiGDh1K48aNARgwYIDDyY2TrCgYUwccP36c4uJimjZtyieffMINN9zArl27AAgJCSEyMpIDBw4QFhbG3Xffzb333utwYlNTWVEwphbatWuXtweQkZHBunXr+Pvf/87UqVMJCwvjwgsv9N4SOmTIEBo1auQ9tkGDBg4mNzWdFQVjarhjx46xbt06RIS4uDgOHTpEREQEqkrjxo2Jiori3nvvJTo6GnBPDb148WKHU5vayoqCMTVQWloaK1euJCMjg88//5yioiLGjh3Lu+++S4sWLZg/fz59+/Zl8ODBhISEOB3X1CF295ExDvrxxx9Zs2YNmZmZ7N+/nz//+c8AjBgxgqysLKKjo73zA8XGxtKhQweHE5u6oKK7j1BVv3wB84D9wMYybW2AdOBrz39bl9n2ILAV2AKM9eUzIiMj1Zia5sUXX9SIiAgVEY2IiNAXX3xRVVVLSkq8+8ydO1ejoqI0ODhYAQW0T58+WlxcrKqqe/fu1aKiIkfym7oPWKNn+Lvqz+U45wOXntI2DVihqj2BFZ7XiEg/YBLQ33PMbBGx0TBT66SmppKYmEhOTg6qSk5ODrfccguRkZF07NiRb775BnCPEzRr1owHHniApUuXsm/fPrKzs72DwB07diQ42K7umsDz6+UjEXEBb6rqAM/rLcAIVd0rIp2AD1W1t4g8CKCqT3j2ew94TFUzKnp/u3xkahqXy0VOTs5p7cHBwSQkJPDYY4/hcrkCH8yYMmrSw2sdVHUvgKcwlK6s3RnILLNfnqftNCKSCCSCLcBtap7c3Nxy20+cOMH8+fMDG8aYKvDn5aOzUd6z9OV2YVQ1RVWjVDWqXbt2fo5lzNk50z9U7B8wprYIdFHY57lshOe/+z3teUBYmf26AHsCnM2Yc5acnExoaOhJbaGhoSQnJzuUyJizE+iisBS4yfP9TcAbZdoniUgjEekK9ARWBzibMecsISGBlJQUIiIiEBEiIiJISUkhISHB6WjG+MRvA80isggYAbQF9gHTgSXAK0A4kAtcp6oHPPsnAbcCxcAvVfWdyj7DBpqNMebsOTLQrKo3nGHT6DPsnwxYH9sYYxxUUwaajTHG1ABWFIwxxnhZUTDGGONlRcEYY4xXrZ4lVUTygdPnFDDGGFORCFUt9+nfWl0UjDHGVC+7fGSMMcbLioIxxhgvKwrGGGO8rCgYUw4RSRaRXSJypJL9rhKRRz3fzxeRawOTEEQkVETeEpH/isgmEZlZZtsvROSWQGUxdYcVBWPKELcg4D9AjA+H/BaY7d9UFfqLqvYBhgIXichlnvZ5wL3OxTK1lRUFU+eIyB9FZGqZ14+JyK9FpJmIrBCRdSKyQUQmeLa7RCRbRGYD64AwVc0sXRCqgs/pBRxX1W/L2fZ7T88hSETGef41v0pEnhaRN8vZf4SIrBSRV0TkKxGZKSIJIrLak7X7qceo6lFV/cDzfaEne5fSbcBOEfGlsBnjZUXB1EWLgZ+VeX098CpwDLhaVS8ARgJ/FZHSBZ56AwtVdaiq+vrsy0W4/xCfRET+BLQHbgFCgDnAZao6HKhoZajBwH3AQGAK0EtVY4DngHsqCiIirYDxuNc+L7UGuNjHn8UYwIqCqYNU9XOgvYicLyKDgQJVzcW9wt8fRGQ9sBz3kq8dPIflqGpm+e94Rp2A/FPaHgFaqeod6n4IqA+wXVV3eLYvquD9slR1r6oeB7YByzztGwDXmQ4SkWDP+z6tqtvLbNoPnO/rD2MMBH6NZmMC5d/AtUBH3D0HgATc/1KPVNUiEdkJNPZs+6EKn/Ej0PKUtiwgUkTaeNYKKW+p2TM5Xub7kjKvS4BgEWkArPW0LVXVRz3fpwBfq+qTp7xfY09GY3xmRcHUVYuBZ3Ev8vQTT1tLYL+nIIwEIs7xM7KBG09pexd4D3hLRMYA/wW6iYhLVXdy8mWts6KqJ4AhZdtE5HHcP9fPyzmkF/BJVT/P1E92+cjUSaq6CWgO7C4zYJwKRInIGty9hv+e6XgR+ZOI5AGhIpInIo+Vs9tHwNAy4xKln/0q7oK01NM0FXhXRFbhXoXwoOczokTkuar+jCLSBUgC+gHrROQLESlbHC7CfZnMGJ/Z3EfGnAMReQr4j6qe8Y+viDRT1SOe4jEL96Wev/s511DgflWd4s/PMXWP9RSMOTd/AEIr2ed2EfkC2IT7Us8cv6dyXzZ7JACfY+oY6ykYY4zxsp6CMcYYLysKxhhjvKwoGGOM8bKiYIwxxsuKgjHGGK//D3KpjjSjaxR9AAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig = plt.figure()\n", + "ax = fig.add_subplot(111)\n", + "\n", + "ax.plot(nessy.variables[var_name]['data'][0, :, 5, 5], \n", + " nessy.lev['data'], \n", + " color='black', marker='o', linestyle='--')\n", + "\n", + "ax.xaxis.tick_top()\n", + "ax.set_ylabel(\"Vertical Levels ({0})\".format(nessy.lev['units']))\n", + "ax.set_xlabel(var_name + \" ({0})\".format(nessy.variables[var_name]['units']))\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Extrapolation options" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "# Do extrapolation below and avobe\n", + "# extrapolate = True \n", + "# extrapolate = ('extrapolate', 'extrapolate')\n", + "\n", + "# Do not extrapolate -> below (bottom); avobe (top)\n", + "# extrapolate = False\n", + "# extrapolate = ('bottom', 'top') # Same behaviour as extrapolate = False\n", + "\n", + "# Only avobe extrapolation:\n", + "extrapolate = ('bottom', 'extrapolate') # (False, True)\n", + "# extrapolate = (False, True)\n", + "# extrapolate = (np.nan, 'extrapolate') # (None, True)\n", + "# extrapolate = (None, True)\n", + "# extrapolate = (0, 'extrapolate')\n", + "\n", + "# Only below extrapolation:\n", + "# extrapolate = ('extrapolate', 'top') # (True, False)\n", + "# extrapolate = ('extrapolate', np.nan) # (True, None)\n", + "# extrapolate = ('extrapolate', 0)\n", + "\n", + "# No extrapolation\n", + "# extrapolate = 0 # (0, 0)\n", + "# extrapolate = np.nan # (np.nan, np.nan)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([ 0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 110, 120,\n", + " 130, 140, 150, 160, 170, 180, 190, 200])" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "level_list = np.arange(0, 210, 10)\n", + "if not ascendant:\n", + " level_list = np.flip(level_list)\n", + "level_list" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "interp_linear = nessy.interpolate_vertical(level_list, kind='linear', extrapolate=extrapolate)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "CPU times: user 135 ms, sys: 2.57 ms, total: 138 ms\n", + "Wall time: 137 ms\n" + ] + } + ], + "source": [ + "%time interp_quad = nessy.interpolate_vertical(level_list, kind='quadratic', extrapolate=extrapolate)" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3deXhU9dXA8e/JAiGEHVQgJIOKLIpsEREBQaTuaNRWMbW41LwVsfpWuygqFhvfatVqnxZrrNYtClqJgrugqEhdgqhBdiQJYV8MJATIdt4/7iRkmSSTZCZ3kpzP88yTmd/c5WQI98z9raKqGGOMMQBhbgdgjDEmdFhSMMYYU8GSgjHGmAqWFIwxxlSwpGCMMaaCJQVjjDEV2lxSEJHzRGSdiGwUkT+4HY+bRKSfiHwkImtE5HsRudXtmEKBiISLyEoRedPtWNwmIl1F5D8istb7d3KG2zG5RUT+1/v/ZJWIvCwiUW7HFAxtKimISDjwD+B8YAgwTUSGuBuVq0qA21V1MDAGuLmNfx7lbgXWuB1EiHgceFdVBwHDaKOfi4j0BX4NJKjqKUA4cJW7UQVHm0oKwGhgo6r+oKpFwDzgEpdjco2qblfVr73P83H+w/d1Nyp3iUgscCHwL7djcZuIdAYmAE8DqGqRqua5G5WrIoAOIhIBRAPbXI4nKNpaUugLbKn0Opc2fhEsJyIeYATwhbuRuO4x4HdAmduBhIDjgd3Av73Vaf8SkY5uB+UGVd0KPAzkANuB/ar6vrtRBUdbSwrio6zNz/MhIjHAa8BtqnrA7XjcIiIXAbtUdYXbsYSICGAk8ISqjgAOAm2yHU5EuuHUKvQH+gAdReTn7kYVHG0tKeQC/Sq9jqWV3gL6S0QicRJCmqoucDsel50JTBWRLJyqxbNF5EV3Q3JVLpCrquV3j//BSRJt0TnAZlXdrarFwAJgrMsxBUVbSwpfAQNEpL+ItMNpKFrockyuERHBqS9eo6qPuh2P21T1TlWNVVUPzt/Gh6raKr8N+kNVdwBbRGSgt2gysNrFkNyUA4wRkWjv/5vJtNJG9wi3A2hOqloiIjOB93B6Dzyjqt+7HJabzgSuATJF5Btv2V2q+raLMZnQcguQ5v0S9QNwncvxuEJVvxCR/wBf4/TaWwmkuhtVcIhNnW2MMaZcW6s+MsYYUwdLCsYYYypYUjDGGFPBkoIxxpgKbTYpiEiy2zGEEvs8jrLPoir7PKpq7Z9Hm00KQKv+h20E+zyOss+iKvs8qmrVn0dbTgrGGGOqadHjFHr27Kkej6dR++7evZtevXoFNqAWzD6Po+yzqMo+j6paw+exYsWKParq85do0SOaPR4PGRkZbodhjDEtiohk1/aeVR8ZY4ypYEnBGGNMBUsKxhhjKlhSMMYYU8GSgjHGmApBSwoi0k9EPhKRNSLyvYjc6i3vLiIfiMgG789ulfa5U0Q2isg6ETk3WLEZE2hpmWl4HvMQ9scwPI95SMtMczsk02qlAR6cy7fH+zpwgnmnUALcrqqDgTHAzSIyBGeN1yWqOgBY4n2N972rgJOB84C5IhIexPiMCYi0zDSSFyWTvT8bRcnen03yomRLDCYI0nAGVGfjLC+f7X0duL+1oCUFVd2uql97n+fjLF3XF2fx6+e8mz0HXOp9fgkwT1WPqOpmYCMwOljxGRMody25i8LiwiplhcWFzFoyy6WITOuiwD5gFfAboLDa+4VA4P7WmmXwmoh4gBHAF8CxqrodnMQhIsd4N+sLfF5pt1xvWfVjJeOdeyQuLi54QRtTh9KyUpZvWU762nRy9uf43Ka2cmMcCvwIbAe2eR++nm8HjtRzrMD9rQU9KYhIDPAacJuqHnDWvPa9qY+yGnNwqGoq3rVRExISWu4cHabFOVJyhA83f8iCNQtYuH4huw7uol14OzpEdOBQyaEa28d1sS8tbZMCedR9oa/rYt8F6AP0BsZVet4H+DWw08c+gftbC2pSEJFInISQpqoLvMU7RaS39y6hN7DLW54L9Ku0eyzOJ2eMa/KP5PPOxndIX5vOW+vfIr8on07tOnHBgAu4bPBlnH/i+Sxcv5DkRclVqpCiI6NJmZziYuQm8Mov9vV9s99G7Rf78ov7uErPy3+WP4+uI4ZinIqSylVI0UDg/taClhTEuSV4Glijqo9WemshMB34s/fnG5XKXxKRR3E+nQHAl8GKz5ja7Cncw8J1C0lfm84Hmz7gSOkRekX34sqTryRxcCKT+0+mfUT7iu2ThiYBMGvJLHL25xDXJY6UySkV5SbUKbCf+r/VbwMO+9i/M0cv6GPxfaHvDXQMQKzlf1OzcKqM4nASQuD+1oI2S6qIjAM+BTKBMm/xXTjtCq/g/DY5wE9VdZ93n1nA9Tg9l25T1XfqOkdCQoLahHgmEHL25/D62tdJX5vOJ9mfUKZlxHeJJ3FQIomDEzmz35mEh1lnuJal/GLvT519zeo/52Jf/eLu64IfiIt98xKRFaqa4PO9ljx1tiUF0xRrdq8hfW066WvTydjm/B2d3OvkikQw4rgR1NEGZlyjwAH8+2bv62LfCd8X9+rf7GOC+Uu4qq6k0KKnzjamIVSVjG0ZpK9NZ8GaBazbuw6A0/uezp8n/5nEwYmc1OMkl6NsyxTIx78G2urdMsG5iJdf2EdT+wW/9V7sA8GSgmkV0jLTfNbpl5SV8Gn2pyxYs4DX171O7oFcwiWciZ6J3DL6Fi4ddCl9O9fo+WxqSKPx9djlF3t/Gmh9Xew74vRO783Ri72vb/adGv5rmRqs+si0eOUjiiv3/mkf3p7RfUezevdq9h7aS1REFOeecC6XDb6Mi066iO4dursYcUtTPoq2eo+XVGAq/n2zP+jjuB2p/dt85ed2sQ80qz4yrdqsJbNqjCg+UnqEZTnLSDo1icRBiZx7wrl0bNfyGgTddRBnGoXaRtFeg4+hRDgJo/ybfQK1N9LaxT4UWVIwLdbOgp28se4NsvfXurIgLyS+0IwRtTQFQBbOhT+r0qP89e569lfgL/j+Zm8N9C2VJQXTomz+cXNFQ/HyLctRlIiwCErKSmpsayOKD1D3RX9vte2jgHicmTdHVnr+v/geRRsP3BHQiI37LCmYkKaqZO7KJH2N03X0253fAjDs2GHMPms2iYMTydyZSfKbbXFEcR51X/R/rLZ9B5yLvAenwTa+0msPcAy+v+GXEexRtCZ0WFIwIadMy/g89/OKRLDpx00Iwth+Y3nkJ49w6aBLOb7b8RXbn3rsqSCtbURx+WRpWdR+0d9fbZ+OHL3Aj6XmRb8njavWCf4oWhM6rPeRCQlFpUUszVrKgjULeGPdG+wo2EFkWCSTj59M4qBEpg6cynExx7kdZgApTvVNFrVf+POr7dOJqhf56hf97lhdvvGH9T4yIelg0UHe2/QeC9Ys4M31b7L/yH46Rnbk/AHnkzgokQsHXEiXqC5uh9lIitNQm0XtF/3q3TS74FzcTwAmU/Oi3xW76Jtgs6RgmtW+Q/tYtG4R6WvTeW/TexwuOUz3Dt1JHJzIZYMu45zjz6FDZAe3w/SD4jS+ZlH7Rb/6FAvdcC7uA4GfUPNbf9cgxmuMfywpmICpbVTx1gNbKyabW5q1lFItJbZzLDeOvJHEQYmMjx9PRFhz/Ck2ZFRuGbCD2i/6OdScMbMHzgX+ZOACal70Ozf5NzAm2KxNwQSEr1HFkWGR9Ovcjx/yfgBgUM9BzmRzgxJJ6JPQzJPN+RqV2wGnS+VJ+L7oF1U7Ri9qr8+Px+bUMS2FtSmYoPO1TnFxWTG5+bmknJ1C4qBEBvca3MxRKc4aTquBW6g5KvcQcH+l18fiXOBHAZdR9aIfR0ucItmYhrKkYBqttKyUZTnL6lynuLi0mLvG3xXkSBRn4b7V3seaSj/31bOveLeLw7lzMKZts6RgGuRwyWGW/LCkYp3iPYV7aB/evpnWKS7FqdrxdfEvqLRdD5x6/Z8CQ4DBwHXAVh/HjMNp+DXGQHCX43wGuAjYpaqneMvmc/R/YFcgT1WHi4gH53/2Ou97n6vqr4IVm2mYA0cO8PaGt0lfm87bG96moKiAzu07c+GAC0kclMj5A87njXVvBHCd4iJgI0cv+uUX/rVUXfu2D85F/zrvz/IE0MvHMR/ERuUaU79g3ik8C/wdeL68QFWvLH8uIo9QdUjmJlUdHsR4TAPsOrirYp3ixT8spqi0iGM6HsPVp1xN4uBEJnkmBWCd4kM43wOqX/w34KzIWs6Dc8GfgnPRL7/4N2QMg43KNcYfQe195L0DeLP8TqFSueD8zzxbVTfUtl19rPdRYGXnZVcsT7ksZxllWkb/rv0rlqc8I/aMRq5TnE/Vqp7yBLCZo1MvhwMncvSiX37hH4g18BoTWKHY+2g8sFNVN1Qq6y8iK3GmdrxbVT/1taOIJOPUAxAX19ZnwWwaVWX17tUsWLOA9LXprNyxEoChxwzl7vF3c9ngyzj12FMb0HV0LzW/9a/GaQQu1w7nQn8a8AuOJoATgfYYY9zlVlKYBrxc6fV2IE5V94rIKOB1ETlZVQ9U31FVU3GWfCIhIaHlDrJwSZmW8eXWLysmm9uwz8nLZ8SewUPnPETi4ERO7H5iHUdQnEFd1Rt6V+N0/yzXERgETOLot/4hQH+sf4MxoavZ/3eKSAROJ/BR5WWqegRvC6KqrhCRTTgjiqxuqIF8jSr+2ZCf8XH2x6SvSef1da+zLX8bPx8qfHpde47pCKXal4iwm6lav14GbKFmlc8anCmby3XFudhfTNWLfz8gLOi/rzEmsJq9TUFEzgPuVNWzKpX1AvapaqmIHA98CgxV1To7mVubQlW+RhWHSzjtw9tTWFJIdGQ05514Hv875hjO7PccIpW7kLbHWW83Cufiv5aqE7YdQ9W6/vKfx2GTtBnTsrjSpiAiLwMTgZ4ikgvMVtWngauoWnUEMAGYIyIlOJ3Rf1VfQjA1/WHxH2qMKi7VUkSE9CvT+ckJPyE6Mg8YRs3J2o4Ar+J8wx8M3EjVi3+PYIdvjAkBQUsKqjqtlvJrfZS9BrwWrFhas+352ysmm8s9kFvj/UE9YXzcQS4dlI6zAPvmOo5W3inMGNNWWYtfC7Rx30bS16SzYO0CPs/9HIAB3QfQvUMnTuyez7g4GB8H4+KgZ3T5Xu/idPr6Nc5Arh0+jmy9uYxp6ywptACqyrc7v63oOrpq1yoAxsUN48XE6Uw5oQO9otdRqtlEeNt2N+yFhevgq63tOH/AA0wd+BuO1v33wkb3GmN8saQQokrLSlm+ZXnFYLKsvCyOixF+OXIQz186gSG9fqR9xCrgW5xePiOIkJv4JFu5/b0FrNi+taL30dSB1Uft2uheY4xvtp6CS9LS0pg1axY5OTnExcWRkpLCFVdewYebPyR9bTpvrHudzu13M8kTzk+H9GZ032K6RO307t0BGAOMw6kSGoOzfq8xxtSvrt5HlhRckJaWRnJyMoWFhc4A3xOh3alhnDY5koTYI0z0hDPRE0nXqPKVvXrgJIDyJDACZ0djjGm4UJzmok37/R9/D8MKmXQBjDsDxnvgjNgyYto7M4Cq9kNkPEeTwEBsIJgxpjlYUgi4NAqKbiU6ci85++HR//bg9NjHGdS+B4u//yMdu67ktQ+PMLI3RIZDWRl8txae/Rcs+xTmzduCSKzbv4Qxpo2y6qOASqOk7Hoiwo6u7VtcCrsOQl/vmu2Hi+GrtcKnHymfvgP//S/s904gHh8fT1ZWVvOHbYxpU+qqPrI6iQBSvbNKQgDnbqBbB7h7cQSfr3mYqMjD5Hz3Ail3RvPuu0cTQnR0NCkp1iXUGOMuSwpNVFJSwvLl83hn8ak4E8jVFBUBD3xWypjBtwPtSUpKIjU1lfj4eESE+Ph4UlNTSUqyLqHGGHdZUqhHWloaHo+HsLAwPB4PaWlpHDp0yFvN81/2Fk5h9OnTmDIpk8Ji38fI2V9zreKkpCSysrIoKysjKyvLEoIxJiRYUqhDedfR7OxsVJXs7GxuuOEX3Jjcif1FQ4CxtG+3lL99Ec49H17NgSMPU1JWtavowSL449LIRq5VbIwxzct6H9Vh1qxZzlgCoFs3uPEmmPnrMvodC+v3HuKeD7vTK/p2rh85g65RXb17HVej99E5Jzxez1rFxhgTGqz3US2Ki4tp3749AwYot/4epl8NHaNgyQ/w10Vw3ZRXuGRQIhFhlleNMS2LDV5roKVLP+KFF6/l7aVhnDehlMMlkPYtPP4iZL4C8ZHxXH7rT90O0xhjAs6SQiVbt25k4VvTGDcpg6f/BTsK4N734Z//hN0fAAXerqOp1j5gjGmdgtbQLCLPiMguEVlVqew+EdkqIt94HxdUeu9OEdkoIutE5NxgxeVIo6CoJ2UqZOUJdy/pxutfDaNdxwHclJxBaSeYtaQv7278J/13PkP01/HIQes6aoxp/YLWpiAiE4AC4PnyNZpF5D6gQFUfrrbtEJwlOkcDfYDFwEmqWlrXORrXppDGkeJraR9ZUlGiCgosWgefbBrJhUMeYpLnbERs7WFjTOvjSpuCqn4iIh4/N78EmKeqR4DNIrIRJ0H8N9BxFRTcSkxMSZUyEdh2AGa83Yetv1kR6FMaY0yL4cY4hZki8p23eqmbt6wvVYcD53rLahCRZBHJEJGM3bt3N/jk0dF7fZb37uSsd2yMMW1ZcyeFJ4ATgOHAduARb7mvehqf9VqqmqqqCaqa0KtXrwYHkFPLuvS+Rh0bY0xb06xJQVV3qmqpqpYBT+FUEYFzZ9Cv0qaxwLZgxPDooz04eKhm+Sur4JKBlwTjlMYY02I0a1IQkd6VXiYC5T2TFgJXiUh7EekPDAC+DEYMp5/+ODNnRJK1DcoUsvMg50eYPiyMtzY8z9YDW4NxWmOMaRGC2SX1ZZyG4oEikisiNwAPiUimiHwHTAL+F0BVvwdeAVYD7wI319fzqLGSkpI455x/M3FsPBHhwlnD48lcnsIxMeE8OCWf6964ljItC8apjTEm5NXbJVVEzgB+jrMuZG/gEM43/LeAF1V1f7CDrE0gprkoKCjg6quvZsyYpdx1Vz7TXoPIHdfw/M3PByhKY4wJLY1eZEdE3gF+CbwHnIeTFIYAdwNRwBsiMjWw4Tavl156iUWLFnHPPfksXw5zz4MlB1/gwWcedDs0Y4xpdnXeKYhIT1XdU+cB/NgmWAJxp+DxeMjOzgbgxBMhMxM0HNpHUGWNZZvl1BjTWjT6TqH6xV5EOotI9/KHr21ampxKfVRPOw0kDDpEQpiApyv83+S9LN50HWmZaS5GaYwxzcOvhmYR+R8R2Ql8B6zwPoIzZ3Uzi4s7OjbhgQegfdU1cujYDlImF3P7u7+hJU8zbowx/vC399EdwMmq6lHV/t7H8cEMrLmkpKQQHR0NQFwtY9f6dIJlN+zisU+imL3wfDJ3fNeMERpjTPPxd+6jTUBhMANxS/mMp7NmzSInJxuPp+Y2ewrhh73CzDOLiIx4l5z97/LCt90pKrmYcfG/Y2DPIc0btDHGBIm/dwp3AstF5EkR+Vv5I5iBNaekpCSysrLweF70ucbyb9+PZPfhFygp3krmhjsoLRvKz07+kRtGPUeXqJN5cWVn5qRPIH5EH8LCwvB4PKSlWRuEMabl8WvqbBH5ElgGZAIVI7tU9bnghVa/4CzHmVZjjWVfvY8KC3cyf+F0uvZewpTTS4iJgr2F8MZKWPAGfDq/PXMfeNrWXjDGhJy6eh/5mxSWq+rYgEfWRMFco9lfJSUlHH98b0aO38Pl18DFE6BrNBw4Am99J0TH/JxRfe4mtstJrsZpjDHlApEUUoBsYBFwpLxcVfcFKsjGCIWkABAWFlbRMykyEiZfApf/Ai6dBD1joLAYPv2hI4dLLmR03Bx6dxrocsTGmLYsEElhs49idbsHUqgkhcoD4Co79tge/M9vB9N34OdcNK6EPl3hSAms2N6dAwfPYUTf2Rzb6Wgj9bKcGXi6ptKnUynb8sPJyktmXNzc5vxVjDFtQJOTQqgKlaSQlpZGcnIyhYVHO2hFR0dXrOdcVlZGRsaXrN/5HCeO+JZ+nTPo27mYklJYvrk9O/eeSa9ePTgt9lU6VmrnPlgEK3fcZInBGBNQdSUFVLXWBzCunvc7A6fUtU0wH6NGjdJQ8eKLL2p8fLyKiMbHx+uLL75Y67ZlZaX6xPwbdO5b3XTtDufXKSvz/Wtu2R/eoGMbY0x9gAyt5bpa39xHfwVOx5nOegWwG2civBNxpr6OB25X1a8ClsIaIFTuFJqioCCf15elkHTug4iP9efKFDrFRNd6F2KMMQ3VpOoj7zrKVwBncnTq7DXAW6q6LMCxNkhrSArlcg9EENu55hISpWXw2Dx45hFY/fXR8vj4eLKyspovQGNMq9HoCfEAVPVHVX1KVa9V1XNV9VJVvbO+hCAiz4jILhFZVansLyKyVkS+E5F0EenqLfeIyCER+cb7+GdDf8mWLisvmYNFVcsOl0BGDtxyJXy/Aj5fD8l3QecuVSfyM8aYQAnmcpzP4qzBUNkHOG0QpwLrcUZKl9ukqsO9j18FMa6QNC5uLit33ETugXDKFHIPhJOx7SaunBhP31PhtichuhM8mQI7dsHz7yq3/P4UXnrpRQ4fPux2+MaYViKovY9ExAO8qaqn+HgvEbhCVZPq2q4uran6qDZVejZFwKgr4IZkmHYGdI2Czbsg7bV2dOc2ZtxkCwMZY+rXpOqjILoeeKfS6/4islJEPhaR8bXtJCLJIpIhIhm7d+8OfpQuS0pKIjU1lfj4eKRU2PPfeDpve5HdB7/lpe8uIutgBHffVMSvfvUQG/edwNYfH+Ff//o7+/e7tkqqMaYlq61bUuUH8FOgk/f53cACYKQf+3mAVT7KZwHpHL1TaQ/08D4fBWwBOtd3/FDqkuqWIyVH9O31c/X5b0/UzT86H82+AvSJeWE6697z9JNPPtaysjJVbVi3WWNM60Vju6SWE5HvVPVUERkH/B/wMHCXqp5ez34eqlULich04FfAZFX1OR23iCwF7lDVOuuG2kL1UUNk523mw8330qndq1x40hE6RMK3m+D19O7E97yTm2+eTeEJhUz7LTxwMcR1gR8LYujR6Z+AdW81pq2oq/rI3/UUyvtKXgg8oapviMh9jQjkPOD3wFmVE4KI9AL2qWqpiBwPDAB+aOjx27r4rv25bsQLlJY9y5LNr7Fx759I6JvJ7Dv2UVTyW2JOhB/K4OYxVIyc7tGpgJKy64kIA0sMxhh/7xTeBLYC5+BU7xwCvlTVYXXs8zIwEegJ7ARm4/Q2ag/s9W72uar+SkQuB+YAJTgJaLaqLqovLrtTqN+Ogh28tf7P7N//ONeMhl4da9syHshqvsCMMa4JxIR40TjdSzNVdYOI9AaGqur7gQ21YSwp+C/eE8/2G3I4fDeE+Rg5DUKlpTKMMa1Yo3sfiUh3EemOM7XFUmCv9/URwK7GLcgDKQ9Q8qOQU2unpFoWqDbGtCn1tSmsABTna2R1Crg6dbbxX1JSEp8d+Iy7ljzBC4kQXunrQElZOyLCUtwLzhgTMuq8U1DV/qp6vPdn9YclhBbmnp/fw6urobjMWRmuTGHLjjCmX1OExzPL1pU2xvg3eE0cPxeRe7yv40RkdHBDM4G2aP0iEvpAVAR888Wv6RQTTVzvMl56CbKzs0lOTrbEYEwb5++I5rnAGcDV3tf5wD+CEpEJmvS16Vw2qDsAt922oMp03ACFhYXMmjXLjdCMMSHC36RwuqreDBwGZ+ZUoF3du5hQcuDIAZb8sISLB3YGhvDNN1t9bmezrxrTtvmbFIpFJByncbl8sJn1X2wh0jLT+O37fdlwSzEDe2RxqGg906b57oocF2e9kIxpy/wd0fw3nLmKjhGRFJxFd+4OWlQmYNIy01i86Tr+fkFxxSjmDu1KeOppiIiI4PnnSyq2jY6OJiXFeiEZ05b5PXW2iAwCJuN0T12iqmuCGZg/bPBa/eIfi+fja3PwdK35XkFBD045JYacnBzi4uJISUmxJT6NaQOaPPeRiDwOzFdVa1wOcWlpacyaNYvsfdkcd04XzrxqP/FdfG8bE7OPrKw9zRugMSak+Vt99DVwt4ichFONNL++GUxN83vmhWeY+fdfMWFGMX/8CSQO3k/n9lBSBhE+p7aw9gNjTFV+JQVVfQ54zjvFxeXAgyISp6oDghqdqVeZlrEs5xM+y3mYLgPeIutDOKYj5BXAq6/Cyx9D75/APy85OjMq2ChmY4xv/t4plDsRGISzeM7qgEdj/LZp30be2/gXROZxwYADTBgPh47AosXw0r/gnbehqMi78RfQLaYHvzljL3FdoLC4BzHtHsemyjbGVOdvm8KDwGXAJmA+cL+q5gUzMFNT3uE83t3wD/KLUhnbL4cZo6GkFN5fDikvdWBhurBzZ811i+IPxPO387MqXsfYCBNjTC38vVPYDJyhqtYq2cxKykr4OGsem/f+mVN7f89VQ53ynP0e7v/LEfJ3n8uUKVfzj79N5Kxxr5CcnFxlpLJ1MzXGNIS/6ymE4UxxcbyqzhGROOA4Vf0y2AHWpTV3Sf1+1zIycu4kvsd/meApJUxg5Sb45NO+/Hr6Z4jE+9yvvPeRdTM1xtQmEIvsPIEzgvlsVR0sIt2A91X1tDr2eQa4CNhVvkazt6F6Pk6bRBbwM++UGYjIncANOCuv/VpV36svrpacFJ754BwmJyyhX1fYkgdLMiZzwZhUlqy5jb7HfMzYfgdoFw4btsM7H3VFfryaSWfdxMknn4yIz65Exhjjl0Akha9VdaSIrFTVEd6yb+tZjnMCUAA8XykpPISzFvOfReQPQDdV/b2IDAFeBkYDfYDFwEmqWlrL4YGWmxSe+eAcrjxrSbXeQE77QFQk7CyIYHvBeCILphMdPp7+/W2WcmNM4DR65bVKGjz3kap+AuyrVnwJ8Jz3+XPApZXK56nqEVXdDGzESRCt0uSEqgkBICLMWefguden0654O8OP+5CTT5xuCcEY06z8TQrV5z5aBvxfI853rKpuB/D+PMZb3hfYUmm7XG9ZDSKSLCIZIpKxe/fuRoTgvn4+ppwAZxzB9EufpVu3ns0bkDHGePk7eC1NRFZwdI2Re8sAABdQSURBVO6jSwM891Fty336iiUVSAWn+iiAMTSbLXkQ383/cmOMaS7+3imgqmtV9R+q+ndVXSMijZl4f6eI9Abw/tzlLc8F+lXaLhbY1ojjtwhLMiZzsKhq2cEip9wYY9zkd1LwoTFdYBYC073PpwNvVCq/SkTai0h/YADganfXYLp+ymLmfzyZkjJQhewfYf7Hk7l+ymK3QzPGtHENneaisjqrbkTkZWAi0FNEcoHZwJ+BV0TkBiAH+CmAqn4vIq/gTJ1RAtxcX8+jlu76KYvZnh/Oxn0nMj5+HddPcTsiY4ypJymIyG9qewuIqWtfVZ1Wy1s+60hUNQWwobfGGOOi+u4UOtXx3uOBDMQYY4z76kwKqvrH5gqkrVmWM4MxsWUcF7Oe3AMRZOUlMy5urtthGWPauKY0NJtGWpYzgxHHPUFEGIhAbOdSRhz3BMtyZrgdmjGmjbOk4AJP19QaI5o7tnPKjTHGTZYUXNCnk++OVbWVG2NMc2ls7yMAVPXRwIbTNmzLDye2c80E4JS7EJAxxnjVd6fQqZ6HaYSsvGSfI5qz8pLdCcgYY7ys95ELxsXNZVkOjIl9gnCBrfnh1vvIGBMS/F2jOQpnAZyTgajyclW9PkhxtXrj4uayPf/JihHNVm1kjAkF/jY0vwAcB5wLfIwzYV1+sIIyxhjjDn+Twomqeg9wUFWfAy4EhgYvLGOMMW7we+U17888ETkF6IKzzrJppGU5M+jVsYxxcc6IZhu4ZowJBf4mhVQR6QbcjTPN9WrgoaBF1crZiGZjTKgS1Ra5eBngrLyWkZHhdhgNlnsgwuc4hdwD4cR2LnEhImNMWyIiK1Q1wdd7ft0piMgDItK10utuIvKnQAXY1tiIZmNMqPK3+uh8Vc0rf6GqPwIXBCek1m9bfniDyo0xprn4mxTCRaR9+QsR6QC0r2P7WonIQBH5ptLjgIjcJiL3icjWSuWtNunYiGZjTKjydznOF4ElIvJvnGU4rweea8wJVXUdMBxARMKBrUA6cB3wV1V9uDHHbUlsRLMxJlT5lRRU9SERycRZSlOA+1X1vQCcfzKwSVWzRSQAh2s5bESzMSYU+XungKq+A7wT4PNfBbxc6fVMEfkFkAHc7m27qEJEkoFkgLi4uACHY4wxbVudbQoissz7M99b91/+yBeRA005sYi0A6YCr3qLngBOwKla2g484ms/VU1V1QRVTejVq1dTQjDGGFNNnUlBVcd5f3ZS1c6VHp1UtakVHucDX6vqTu85dqpqqaqWAU8Bo5t4/JBmI5qNMaHI33EKL/hT1kDTqFR1JCK9K72XCKxq4vFDlo1oNsaEKn+7pJ5c+YWIRACjGntSEYkGpgALKhU/JCKZIvIdMAn438YeP9TZGs3GmFBV33KcdwJ3AR0qtSEIUAQ0+gqmqoVAj2pl1zT2eC2NjWg2xoSq+toU/g9nRtTnq7Un9FDVO5snxNbHRjQbY0JVvdVH3obfYc0QS5thI5qNMaHK3zaFz0XktKBG0oaMi5vLyh03UVIGqs7sqCt33GQjmo0xrvNr6mwRWQ0MBLKAgzjtCqqqpwY1unq01Kmzy23PD68Y0WyMMc2lrqmz/R3RfH4A4zHGGBOi/Ko+UtVsoB9wtvd5ob/7GmOMaTn8Hbw2G/g9UN7jKBJn5lTTSDai2RgTivz9tp+IM0/RQQBV3QZ0ClZQrZ2NaDbGhCp/k0KROi3SCiAiHYMXUutnI5qNMaHK36Twiog8CXQVkRuBxTiT1plGsBHNxphQ5e8iOw+LyBTgAE7X1HtV9YOgRtaKbcsPJ7ZzzQTglLsQkDHGeNW3nsLfRWQsgKp+oKq/VdU7LCE0jY1oNsaEqvqqjzYAj4hIlog8KCLDmyOo1s5GNBtjQpW/I5rjcZbOvAqIwlkHYZ6qrg9ueHWzEc3GGNNwdY1o9nvwmqo+qKojgKtxuqiuCWCMxhhjQoC/g9ciReRiEUkD3gHWA5c39qTe6qhMEflGRDK8Zd1F5AMR2eD92a2xx28JbPCaMSYU1dfQPEVEngFygWTgbeAEVb1SVV9v4rknqerwSrcwfwCWqOoAYIn3datkg9eMMaGqvjuFu4D/AoNV9WJVTVPVg0GK5RLgOe/z54BLg3Qe19ngNWNMqKpznIKqTgrSeRV4X0QUeFJVU4FjVXW797zbReQYXzuKSDLOXQtxcXFBCi+4bPCaMSZUuTXT6ZmqOhJnSu6bRWSCvzuqaqqqJqhqQq9evYIXYRDZcpzGmFDlSlLwTqiHqu4C0oHRwE4R6Q3g/bnLjdiagw1eM8aEqmZPCiLSUUQ6lT8HfgKsAhYC072bTQfeaO7YmosNXjPGhCp/V14LpGOBdBEpP/9LqvquiHyFM/HeDUAO8FMXYms24+Lmsj3/yYrBazbnkTEmFDR7UlDVH4BhPsr3ApObOx5jjDFH2ZKaxhhjKrhRfWRwBrCNiS3juBhnRHNWXrK1KRgDFBcXk5uby+HDh90OpcWLiooiNjaWyMhIv/expOCCyiOawRnR3C3qCZblYInBtHm5ubl06tQJj8eDt+3RNIKqsnfvXnJzc+nfv7/f+1n1kQtsRLMxtTt8+DA9evSwhNBEIkKPHj0afMdlScEFNqLZmLpZQgiMxnyOlhRcYCOajTGhypKCC2xEszGBk5aWhsfjISwsDI/HQ1paWrOd+4ILLiAvL6/Obe69914WL17cqOMvXbqUiy66qFH7NpY1NLtgXNxcluXAmNgnCBfYmh9uvY+MaYS0tDSSk5MpLCwEIDs7m+Rk58tVUlJS0M6rqqgqb7/9dr3bzpkzJ2hxBIPdKbhkXNxcdh8MY1nOScR2LrGEYEwtJk6cWOMxd67z/+XOO++sSAjlCgsLufXWWwHYs2dPjX399eijj3LKKadwyimn8Nhjj5GVlcXgwYOZMWMGI0eOZMuWLXg8Hvbs2QPA/fffz6BBg5gyZQrTpk3j4YcfBuDaa6/lP//5DwAej4fZs2czcuRIhg4dytq1awH48ssvGTt2LCNGjGDs2LGsW+feEr2WFIwxLVZubq7P8r179zbpuCtWrODf//43X3zxBZ9//jlPPfUUP/74I+vWreMXv/gFK1euJD4+vmL7jIwMXnvtNVauXMmCBQuoa+34nj178vXXX3PTTTdVJI5BgwbxySefsHLlSubMmcNdd93VpPibwqqPjDEhbenSpbW+FxcXR3Z2do3y8gt2z54969y/NsuWLSMxMZGOHTsCcNlll/Hpp58SHx/PmDFjfG5/ySWX0KFDBwAuvvjiWo992WWXATBq1CgWLFgAwP79+5k+fTobNmxARCguLm5wzIFidwousTWajWm6lJQUoqOjq5RFR0eTkpLSpOOqqs/y8iTh7/a+tG/fHoDw8HBKSkoAuOeee5g0aRKrVq1i0aJFro7mtqTgAluj2ZjASEpKIjU1lfj4eESE+Ph4UlNTm9zIPGHCBF5//XUKCws5ePAg6enpjB8/vtbtx40bV3ExLygo4K233mrQ+fbv30/fvn0BePbZZ5sSepNZ9ZEL6h7RbA3OxjREUlJSwHsajRw5kmuvvZbRo0cD8Mtf/pJu3brVuv1pp53G1KlTGTZsGPHx8SQkJNClSxe/z/e73/2O6dOn8+ijj3L22Wc3Of6mkIbc9oSahIQEratBJ1SVqRDmY6BhmUKYtNx/D2MCYc2aNQwePNjtMBqsoKCAmJgYCgsLmTBhAqmpqYwcOdLtsHx+niKyQlUTfG1vdwou2JYfTmznmlNaOOUuBGSMabLk5GRWr17N4cOHmT59ekgkhMZo9qQgIv2A54HjgDIgVVUfF5H7gBuB3d5N71LV+keGtEBZecl0i3qiShVS+YhmSwrGtEwvvfSS2yEEhBt3CiXA7ar6tXet5hUi8oH3vb+q6sMuxNSsbESzMSZUubEc53Zgu/d5voisAfo2dxxuszWajTGhyNUuqSLiAUYAX3iLZorIdyLyjIj4bOoXkWQRyRCRjN27d/vaxBhjTCO5lhREJAZ4DbhNVQ8ATwAnAMNx7iQe8bWfqqaqaoKqJvTq1avZ4jXGmLbAlaQgIpE4CSFNVRcAqOpOVS1V1TLgKWC0G7E1FxvRbExgpGWm4XnMQ9gfw/A85iEts+lTZ8fExACwbds2rrjiiiYfryVp9qQgzlJATwNrVPXRSuW9K22WCKxq7tiai41oNiYw0jLTSF6UTPb+bBQle382yYuSA5IYAPr06VMxw2mwlE91ESrc6H10JnANkCki33jL7gKmichwQIEs4H9ciK1Z2IhmY/xz27u38c2Ob2p9//PczzlSeqRKWWFxITe8cQNPrXjK5z7DjxvOY+c95tf5s7KyuOiii1i1ahXPPvssCxcupLCwkE2bNpGYmMhDDz0EwPvvv8/s2bM5cuQIJ5xwAv/+97+JiYlhzpw5LFq0iEOHDjF27FiefPJJRISJEycyduxYPvvsM6ZOncrtt9/u5ycSfM1+p6Cqy1RVVPVUVR3ufbytqteo6lBv+VRvL6VWydZoNiYwqieE+sqb6ptvvmH+/PlkZmYyf/58tmzZwp49e/jTn/7E4sWL+frrr0lISODRR51KkJkzZ/LVV1+xatUqDh06xJtvvllxrLy8PD7++OOQSghgI5pdYSOajfFPfd/oPY95yN7vY+rsLvEsvXZpwOOZPHlyxZxGQ4YMITs7m7y8PFavXs2ZZ54JQFFREWeccQYAH330EQ899BCFhYXs27ePk08+uWJa7SuvvDLg8QWCJQUX2IhmYwIjZXIKyYuSKSw+uvpadGQ0KZObNnV2bcqnvYajU1+rKlOmTOHll1+usu3hw4eZMWMGGRkZ9OvXj/vuu6/KlNi1TcPtNps62wXj4uaycsdNlJSBKuQeCGfljptsRLMxDZQ0NInUi1OJ7xKPIMR3iSf14lSShgZvfebqxowZw2effcbGjRsBZznQ9evXVySAnj17UlBQEPQG60CxOwWX2IhmYwIjaWhSsyaB6nr16sWzzz7LtGnTOHLEacv405/+xEknncSNN97I0KFD8Xg8nHbaaa7F2BA2dbaLtueHVyQFY4yjpU6dHaoaOnW2VR8ZY4ypYEnBJTai2RgTiiwpuMBGNBtjQpUlBRfUPaLZGGPcY0nBBTai2RgTqiwpuGBbfniDyo0xprlYUnBBVl4yB4uqlpWPaDbGNFQa4MG5nHm8r0PTtdde2+BBbK+//jqrV6+ueH3vvfeyePHiQIdWwZKCC2xEszGBkgYkA9k4Eyxne1+HbmLwpbS09qrj6klhzpw5nHPOOUGLxZKCS8bFzWX3wTCW5ZxEbOcSSwjG+HQbMLGOxw1AYbV9Cr3lte1zm19nTklJYeDAgZxzzjlMmzaNhx9+mIkTJ1I+YHbPnj14PB7AmWJ7/PjxjBw5kpEjR7J8+XIAVJWZM2cyZMgQLrzwQnbt2lVxfI/Hw5w5cxg3bhyvvvoqTz31FKeddhrDhg3j8ssvp7CwkOXLl7Nw4UJ++9vfMnz4cDZt2lTlbuOrr75i7NixDBs2jNGjR5Ofn+/X71YXm+bCGNOC1TZFdtOmzl6xYgXz5s1j5cqVlJSUMHLkSEaNGlXr9scccwwffPABUVFRbNiwgWnTppGRkUF6ejrr1q0jMzOTnTt3MmTIEK6//vqK/aKioli2bBkAe/fu5cYbbwTg7rvv5umnn+aWW25h6tSpXHTRRTVWgCsqKuLKK69k/vz5nHbaaRw4cIAOHTo06fcGSwrGmJBW32I4Hpwqo+rigaWNPuunn35KYmIi0dHRAEydOrXO7YuLi5k5cybffPMN4eHhrF+/HoBPPvmEadOmER4eTp8+fTj77LOr7Fd5+uxVq1Zx9913k5eXR0FBAeeee26d51y3bh29e/eumFOpc+fATKAWctVHInKeiKwTkY0i8ge34wkWG9FsTCCkANHVyqK95U3jrBxcVUREBGVlZQBVpsH+61//yrHHHsu3335LRkYGRUVFdR6nXOXps6+99lr+/ve/k5mZyezZs6sc3xdVrfPYjRVSSUFEwoF/AOcDQ3CW6BziblSBZyOajQmUJCAV585AvD9TveWNN2HCBNLT0zl06BD5+fksWrQIcNoBVqxYAVClF9H+/fvp3bs3YWFhvPDCCxUNxxMmTGDevHmUlpayfft2Pvroo1rPmZ+fT+/evSkuLiYt7WhDeadOnXy2FQwaNIht27bx1VdfVewfiPWeQyopAKOBjar6g6oWAfOAS1yOKeBsRLMxgZSEs6x7mfdn06fRHjlyJFdeeSXDhw/n8ssvZ/z48QDccccdPPHEE4wdO5Y9e/ZUbD9jxgyee+45xowZw/r16yvuABITExkwYABDhw7lpptu4qyzzqr1nPfffz+nn346U6ZMYdCgQRXlV111FX/5y18YMWIEmzZtqihv164d8+fP55ZbbmHYsGFMmTKl3rsLf4TU1NkicgVwnqr+0vv6GuB0VZ1ZaZtknD5nxMXFjcrO9lWfGNrKVAjzcddXphAmofPvYYwbQnHq7Pvuu4+YmBjuuOMOt0NpsJY+dbavCrIqV0lVTVXVBFVN6NWrVzOFFVg2otkYE6pCrfdRLtCv0utYYJtLsQSNrdFsTMty3333uR1Cswm1O4WvgAEi0l9E2gFXAQtdjingykc05x4Ip8xGNBtTQyhVa7dkjfkcQ+pOQVVLRGQm8B4QDjyjqt+7HFZQOAnASQKxnbE7BGO8oqKi2Lt3Lz169AhKl8u2QlXZu3cvUVFRDdovpJICgKq+DbztdhzGGHfExsaSm5vL7t273Q6lxYuKiiI2NrZB+4RcUjDGtG2RkZH079/f7TDarFBrUzDGGOMiSwrGGGMqWFIwxhhTIaRGNDeUiOzG9xSJxhhjahevqj5H/7bopGCMMSawrPrIGGNMBUsKxhhjKlhSMMYYU8GSgjE+iEiKiGwRkYJ6trtURO71Pn/WO/17sxCRaBF5S0TWisj3IvLnSu/NFJHrmisW03pYUjCmEnGEAYtwFn2qz+8on8TKHQ+r6iBgBHCmiJzvLX8G+LV7YZmWypKCaXVE5EERmVHp9X0icruIxIjIEhH5WkQyReQS7/seEVkjInOBr4F+qvq5qm6v5zwnAUdUdY+P9+733jmEicgF3m/zy0TkbyLypo/tJ4rIxyLyioisF5E/i0iSiHzpjfWE6vuoaqGqfuR9XuSNPbb8PSBLRPxJbMZUsKRgWqN5wJWVXv8MeBU4DCSq6khgEvCIHJ2GcyDwvKqOUFV/x76ciXMhrkJEHgKOAa4D2gFPAuer6jigrpWhhgG3AkOBa4CTVHU08C/glroCEZGuwMXAkkrFGcB4P38XYwBLCqYVUtWVwDEi0kdEhgE/qmoOzsp+D4jId8BioC9wrHe3bFX9vIGn6g1Un8rzHqCrqv6POoOABgE/qOpm7/sv13G8r1R1u6oeATYB73vLMwFPbTuJSIT3uH9T1R8qvbUL6OPvL2MM2CyppvX6D3AFcBzOnQM4K7r3AkaparGIZAHlk80fbMQ5DgFdqpV9BYwSke6qug/fS8zW5kil52WVXpcBESISDqzwli1U1Xu9z1OBDar6WLXjRXljNMZvlhRMazUPeAroCZzlLesC7PImhElAfBPPsQb4ebWyd3EWiXpLRH4CrAWOFxGPqmZRtVqrQVS1FBheuUxE/oTze/3Sxy4nAZ819nymbbLqI9MqeVfs6wRsrdRgnAYkiEgGzl3D2tr2F5GHRCQXiBaRXBG5z8dmnwAjKrVLlJ/7VZyEVL6U7AzgXRFZBuwE9nvPkSAi/2rs7ygiscAsYAjwtYh8IyKVk8OZONVkxvjN5j4ypglE5HFgkarWevEVkRhVLfAmj3/gVPX8NchxjQB+o6rXBPM8pvWxOwVjmuYBILqebW4UkW+A73Gqep4MelROtdk9zXAe08rYnYIxxpgKdqdgjDGmgiUFY4wxFSwpGGOMqWBJwRhjTAVLCsYYYyr8P0U38EXFQQwHAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig = plt.figure()\n", + "ax = fig.add_subplot(111)\n", + "\n", + "ax.plot(nessy.variables[var_name]['data'][0, :, 5, 5], \n", + " nessy.lev['data'], \n", + " color='black', marker='o', linestyle='--', label='original')\n", + "\n", + "ax.plot(interp_linear.variables[var_name]['data'][0, :, 5, 5], \n", + " level_list, \n", + " color='green', marker='o', label = 'linear')\n", + "\n", + "ax.plot(interp_quad.variables[var_name]['data'][0, :, 5, 5], \n", + " level_list, \n", + " color='yellow', marker='o', label = 'quadratic')\n", + "\n", + "ax.xaxis.tick_top()\n", + "ax.set_ylabel(\"Vertical Levels ({0})\".format(nessy.lev['units']))\n", + "ax.set_xlabel(var_name + \" ({0})\".format(nessy.variables[var_name]['units']))\n", + "\n", + "ax.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.7.4" + }, + "toc-autonumbering": true, + "vscode": { + "interpreter": { + "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6" + } + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/tutorials/6.Others/6.5.Add_new_variables.ipynb b/tutorials/6.Others/6.5.Add_new_variables.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..04c5c487803267fc48de6016d5806b515f35428b --- /dev/null +++ b/tutorials/6.Others/6.5.Add_new_variables.ipynb @@ -0,0 +1,192 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# How to add new variables" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "from nes import *\n", + "import numpy as np\n", + "import matplotlib\n", + "import matplotlib.pyplot as plt\n", + "import geopandas as gpd" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "original_file_path = \"/gpfs/projects/bsc32/models/NES_tutorial_data/MONARCH_d01_2022111512.nc\"" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Open original file" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "# Open file with only the metadata\n", + "nessy = open_netcdf(original_file_path, parallel_method='Y')" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# Load variables data\n", + "nessy.keep_vars(['NO2']) # If you avoid this part all variables will be loaded\n", + "nessy.load()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Create new variables" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2.1. Create variable my modifying an existen one" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [], + "source": [ + "from copy import deepcopy\n", + "\n", + "# Copy the base variable and all the metadata\n", + "nessy.variables['new_var1'] = deepcopy(nessy.variables['NO2'])\n" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "# Modify the data and the metadata\n", + "nessy.variables['new_var1']['data'] *= 2\n", + "nessy.variables['new_var1']['long_name'] = 'Double NO2'\n", + "nessy.variables['new_var1']['standard_name'] = 'Invented_name'" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "# Create the new file\n", + "nessy.to_netcdf('new_var1.nc')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 2.1. Create variable from scratch" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "nessy = open_netcdf(original_file_path, parallel_method='Y')" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "# Remove all variables information to keep only the metadata\n", + "nessy.variables = {}" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "# Create metadata\n", + "nessy.variables['new_var2'] = {'data': None,\n", + " 'units': 'kg.m-2',\n", + " 'description': 'new variable from scratch',\n", + " 'short_name': 'newvar2',\n", + " }\n", + "\n", + "# Fill the data\n", + "shape = (len(nessy.time), len(nessy.lev['data']), nessy.lat['data'].shape[0], nessy.lon['data'].shape[-1])\n", + "nessy.variables['new_var2']['data'] = np.random.rand(*shape)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "nessy.to_netcdf('new_var2.nc')" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "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.7.4" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}