From d2b45424f61d58c15bc58e9a62a084c574c196ed Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 19 Jan 2023 12:33:40 -0600 Subject: [PATCH 01/14] Changes to optionally return distances with .sel Changes include: * accessor optionally takes in a "distances_name" str with which to name the new variable in results to store calculated distances from query() * query() now returns distances * distances are returned (always?) in radians from the trees, but are converted to kilometers in query. * If a DataArray was input to .sel with "distances_name" input, a Dataset will be returned to contain the new variable * added example to intro notebook * added a test to test_accessor * also fixed an import in conftest Notes: * I did not test this with the dask approach --- doc/examples/introduction.ipynb | 234 ++++++++++++++++++++++++++++---- src/xoak/accessor.py | 37 ++++- src/xoak/tests/conftest.py | 2 +- src/xoak/tests/test_accessor.py | 21 +++ 4 files changed, 263 insertions(+), 31 deletions(-) diff --git a/doc/examples/introduction.ipynb b/doc/examples/introduction.ipynb index 082a960..6a1b82f 100644 --- a/doc/examples/introduction.ipynb +++ b/doc/examples/introduction.ipynb @@ -11,7 +11,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "metadata": {}, "outputs": [], "source": [ @@ -31,7 +31,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 2, "metadata": {}, "outputs": [], "source": [ @@ -44,9 +44,37 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 3, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/html": [ + "
<xarray.Dataset>\n",
+       "Dimensions:  (x: 100, y: 100)\n",
+       "Coordinates:\n",
+       "    lat      (x, y) float64 -46.35 76.28 -85.01 -13.43 ... 62.3 -88.36 5.12\n",
+       "    lon      (x, y) float64 -1.377 53.78 52.71 -174.7 ... -164.6 -116.2 -51.76\n",
+       "Dimensions without coordinates: x, y\n",
+       "Data variables:\n",
+       "    field    (x, y) float64 -47.72 130.1 -32.3 -188.1 ... -102.3 -204.5 -46.64
" + ], + "text/plain": [ + "\n", + "Dimensions: (x: 100, y: 100)\n", + "Coordinates:\n", + " lat (x, y) float64 -46.35 76.28 -85.01 -13.43 ... 62.3 -88.36 5.12\n", + " lon (x, y) float64 -1.377 53.78 52.71 -174.7 ... -164.6 -116.2 -51.76\n", + "Dimensions without coordinates: x, y\n", + "Data variables:\n", + " field (x, y) float64 -47.72 130.1 -32.3 -188.1 ... -102.3 -204.5 -46.64" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "ds_mesh = xr.Dataset(\n", " coords={'lat': (('x', 'y'), lat), 'lon': (('x', 'y'), lon)},\n", @@ -67,7 +95,7 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 4, "metadata": {}, "outputs": [], "source": [ @@ -83,9 +111,33 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 5, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/html": [ + "
<xarray.Dataset>\n",
+       "Dimensions:    (trajectory: 30)\n",
+       "Dimensions without coordinates: trajectory\n",
+       "Data variables:\n",
+       "    latitude   (trajectory) float64 -10.0 -8.276 -6.552 ... 36.55 38.28 40.0\n",
+       "    longitude  (trajectory) float64 -150.0 -139.7 -129.3 ... 129.3 139.7 150.0
" + ], + "text/plain": [ + "\n", + "Dimensions: (trajectory: 30)\n", + "Dimensions without coordinates: trajectory\n", + "Data variables:\n", + " latitude (trajectory) float64 -10.0 -8.276 -6.552 ... 36.55 38.28 40.0\n", + " longitude (trajectory) float64 -150.0 -139.7 -129.3 ... 129.3 139.7 150.0" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "ds_trajectory = xr.Dataset({\n", " 'latitude': ('trajectory', np.linspace(-10, 40, 30)),\n", @@ -106,9 +158,37 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 6, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/html": [ + "
<xarray.Dataset>\n",
+       "Dimensions:  (trajectory: 30)\n",
+       "Coordinates:\n",
+       "    lat      (trajectory) float64 -11.46 -8.57 -6.349 ... 37.13 38.05 41.56\n",
+       "    lon      (trajectory) float64 -150.6 -139.7 -129.9 ... 129.0 140.2 149.2\n",
+       "Dimensions without coordinates: trajectory\n",
+       "Data variables:\n",
+       "    field    (trajectory) float64 -162.1 -148.3 -136.3 ... 166.2 178.2 190.8
" + ], + "text/plain": [ + "\n", + "Dimensions: (trajectory: 30)\n", + "Coordinates:\n", + " lat (trajectory) float64 -11.46 -8.57 -6.349 ... 37.13 38.05 41.56\n", + " lon (trajectory) float64 -150.6 -139.7 -129.9 ... 129.0 140.2 149.2\n", + "Dimensions without coordinates: trajectory\n", + "Data variables:\n", + " field (trajectory) float64 -162.1 -148.3 -136.3 ... 166.2 178.2 190.8" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "ds_selection = ds_mesh.xoak.sel(\n", " lat=ds_trajectory.latitude,\n", @@ -127,14 +207,79 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 7, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEGCAYAAABLgMOSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxmUlEQVR4nO3deXhcd3no8e97zizSaLMdSV6jyDZ2VhKHOCYQCikEYgIh7ARyIVRQXyi5lF64kJDrp2lLblNa0lCWUJOqJC0lpLQhYcnOEiCrA47jxHHiRZHlVbK8SBotM+e89485Uka2ZI+k0ZyZ0ft5nvNo5ndmeTUz0ju/XVQVY4wxZqKcsAMwxhhTmiyBGGOMmRRLIMYYYybFEogxxphJsQRijDFmUiJhB5AP9fX12tzcHHYYxpgS8PTTT3epasNUHmNJRZ0m/XROt92XSt6vqqun8nzFqiwSSHNzM+vXrw87DGNMCRCRl6f6GP1+mpaGM3O67d/ufqp+qs9XrMoigRhjTKE5ImGHEDpLIMYYM0ECuJY/LIEYY8xECYJrNRBLIMYYM1EiEHUsgVgCMcaYSbAmLEsgxpgSNNDTy2BfkkMde6lf0kSsqpJoPF6w5xfAwTKIJRBjTEkZ7Euy4a77uWft1/BSaRJz6viT2/+Rk889Eylgv4TNwrbXwBhTYtKDKe694Zt4qcxEvmT3Ye694ZskDx4uWAyZUViS01HOrAZijCkp6aEhBnr6RpUd3rOfQu9tZH0gVgMxxpSYaDxG86oVo8pe/c63EKusKFgMw6OwcjnKmdVAjDElJTG7jiv/+f/xy3/6Hnu3bOO0N7+e8z9yObFEZcFisImEGaEnEBFxgfXALlV9p4jMAX4INANtwAdV9WB4ERpjik3t3HpWf/nPSPb0QtSlLzVIpWphO9EtgYSfQIA/BzYDtcH1a4CHVfVGEbkmuP6lsIIzxhQfVeX5F7fQ2tpKe3s7TU1NtLS0sGLFioIkEZuJnhFqH4iILALeAdyaVXw5cFtw+Tbg3QUOyxhT5Do7O2ltbaWtrQ3f92lra6O1tZXOzs6CPL9Ipgkrl6OchV0DuRn4IlCTVTZXVfcAqOoeEWkc644isgZYA9DU1DTNYRpjikkqlaK9vX1UWXt7O6lUqiDPL0DUaiDh1UBE5J3AflV9ejL3V9V1qrpSVVc2NExpbxhjTImJRqPHfHFsamoiGo0WLAargYTbhHUh8C4RaQPuAN4sIv8O7BOR+QDBz/3hhWiMKUYNDQ20tLTQ3NyM4zg0NzfT0tJCob5M2kTCjNCasFT1WuBaABG5CPiCqv4PEfl74CrgxuDn3WHFaIwpTiLCihUrWLt2LalUimg0SkNDg43CKrBinEh4I/BWEXkJeGtw3RgzwwwcPkLPvk4O795LT9eBY2aaiwiNjY0sXLiQxsbGgiaPTCd6fmogItIqIvtFZFNW2fUisktENgTHpVnnrhWRrSKyRUQumaZfMSdhd6IDoKq/An4VXD4AvCXMeIwx4Ro4fIRffOWfePK2O0mnUpx1+SVc+tX/S03DSQVNFMeTx/6N7wHfBG4/qvwfVfUfsgtE5AzgCuBMYAHwkIgsV1Uvb9FMQDHWQIwxM9zuZ57n8dYfMDQ4hO8rG++6j833/YKenp6wQwMy80AijpPTcSKq+gjQneNTXw7coaqDqroD2AqsmvxvMjWWQIwxoRg80sPg4bETQufmrXje6C/VXZu34nuhfNE+loA4ktMxBVeLyMagiWt2ULYQ2Jl1m46gLBSWQIwxBZVK9rP/98/y6HU38uh1f8u+pzeSSvaPus3SN184aoMoEWHp296U6XwoAgI4ruR0APUisj7rWJPDU9wCLAVWAHuAr2U99dEKuwxxlqLoAzHGzBz9+7u4/+Ofwx8aAqDj14/zrrv/lbolp4zcpmZ+I1d872Z+8/VbSQ0Mcv4nr6D+1CXU1dWFFfZoApJ7J0iXqq6cyMOr6r6RpxL5LvDT4GoHcHLWTRcBuyfy2PlkCcQYU1A7fv4w/tAQvu+jCr4/yPafPsi5n/3kyG3i1VUsf9sbWfCas1BVpCJOTV1t0XSgw4QSyMQfW2T+8IocwHuA4RFa9wD/ISI3kelEXwY8OW2BnIAlEGNMQdWcvJBUKkUymcTzPFzXpXJeYyZRZCUIx3WpnVukq0zIlPs3sh5KfgBcRKapqwP4S+AiEVlBpnmqDfifAKr6nIjcCTwPpIHPhDUCCyyBGGMKbN6FK5mz4kx6fvN45vr5K5jz2hV0dXUVbCb5VImAG3Xz8liq+uExiv/lOLe/AbghL08+RZZAjDEF1ZtOce5X/g/nHOoBVdw5s/jat7/FZz/72bBDm5B81UBKmSUQY0zB3fydWxARRITt27cXfCHEKRNBXBvEaq+AMaagGhoa+PjHP47neSPJo5ALIebDBIfxli2rgRhj8k5V6ezsHHOhw2JYCHHKBMS+flsCMcZMTqqnF3EcIlWJUeWqyoYNG4673ezwQoglSwQnT53opcxyqDFmQtLJJId+/wwv/t1NvHTTP9G3fQfe4ODI+bC3my0ECSYS5nKUM6uBGGMmpL+9g42fvxYN1qU68OgTnP9vt+IGS4+Evd1soTg5LJRY7uwVMMbkzB8aYu/P7x9JHgBeX5LuJ9aPXC+G7WanneRW+7AaiDHGDHNd4mP0XcQb6kcuNzQ08MlPfpJnNmygpjJByvc47YwzSmqU1YkI4Ng8EKuBGGNy57gu897xNqqWLB4pm3PB+dScvnzkuohwxuIlvHV+E0uf3coFboLTT2kurVFWJ2J9IIDVQIwxExSpq+PVN99I8uWduBUVVC6YR7SmZuS8NzBI+w9/zNZvtQKw5ycPcHjDJk770v8iWlsz3sOWFgEnat+/LYEYY0bx02m8ZDIzRLe6etQ5VeWZZ57h9ttvJx6Pc9JJJ3HxxRePGqLrDQzScdfPRt1vz/2/4NQvXl2w32H62Ux0sARijMmS7uvj4KO/o+tXvyRSW8uCD36IykWLcKIxYPQQ3WHt7e2sXbv2lXkdqrgVFaMe14nHkTH3QipNItYHAtYHYowJqO9z+A+/5+V136HvxS0cXv8UL/719XhZuwXmMkQ3Ul3F0jVXIZFXJtot+ZMPj7peDgqwpW3RsxqIMQYAL5nk4KO/G13W20ty+zbqzn0N8MoQ3ewayNFDdJ1ohIY/ei1v+O/b6F6/gdozTqVy4bxjZqyXNBEca8KyBGKMyXBiUeLz5o8uFCGWNWy3oaGBlpaWY5YpOXqIbqQqQaQqQaJpYSFCLzhhenckLBWWQIwxADixOHPfeRk9zz1L79atiOMw953vIjpr9shtymIhxHwQyq5JbjIsgRhjRkRmzeKUL15L7759OJUV7O3u5vD27Zx11lnlsxBiPlgTFmAJxBiTpaurixtuuIGhoSH6+/s5cOAAzc3No0dZGQDE1sKyBGKMeUUqlRpZRXdYOS6EOFViOxICNozXGJNlRiyEmCf5GsYrIq0isl9ENmWVzRGRB0XkpeDn7Kxz14rIVhHZIiKXTNOvl5PQEoiIVIjIkyLyjIg8JyJ/FZSP+8IZY6bX8Cir5uZmHMehubm55LabLQgRnGg0pyMH3wNWH1V2DfCwqi4DHg6uIyJnAFcAZwb3+baIhNabH2YT1iDwZlXtFZEo8FsRuRd4L5kX7kYRuYbMC/elEOM0ZsawUVa5yazGm5/v36r6iIg0H1V8OXBRcPk24Fdk/g9eDtyhqoPADhHZCqwCHstLMBMUWgJRVQV6g6vR4FDGf+GMMXmgqiQPH0I9Dy8Spba2dlSCsFFWORAg93kg9SKyPuv6OlVdd4L7zFXVPQCqukdEht+QhcDjWbfrCMpCEWonelD1ehp4FfAtVX1CRMZ74Y6+7xpgDXBMm60xZmx+Ok3//r3s/cmPGDp8mFlvuAh91enUWS1jYibWid6lqivz9cxjlGmeHnvCQk0gquoBK0RkFnCXiJw1gfuuA9YBrFy5MrQX0JhSkk728cINX2agqwuA7qefYMlffBmJx6mrqws5utIiMq1dyPtEZH7wJXo+sD8o7wBOzrrdImD3dAZyPEUxCktVD5FpqlpN8MIBHPXCGWNyoL6Pl+wlfagbL9mLnxoaOdfftm0keQw7/NgjOGkbpjsxmRpILsck3QNcFVy+Crg7q/wKEYmLyGJgGfDklH6VKQitBiIiDUBKVQ+JSCVwMfB3vPLC3cjoF84Yk4PUgf3s//53SO3fg1NVQ/17/geVy87AicaI1s3Gjbh46Vf2NI/NOQlxbUrYRIiAG83P4CcR+QGZft96EekA/pLM/787ReQTQDvwAQBVfU5E7gSeB9LAZ4KWnFCE+amZD9wW9IM4wJ2q+lMReYwxXjhjzIl5yT66f/6fpPbvAcDv66Hrx//Oos9dD9EY8YZGFrx5Nbt/cR9e2iOxcBHz3345FdZ8NTEieZuJrqofHufUW8a5/Q3ADXl58ikKcxTWRuDcMcoPMM4LZ4w5AfUZ3LcLz/eCoaYufm8POjQECXATVSy44uM0rH4X6d4e4gsWEamqsg70SbCZ6LaUiTFlQ1U5eKQHOXkph3dsw3VdEokEiYWnILH4yO0iVVVEqqpCjLQMiCUQKJJOdGPM1HV2dvKdf2ml8g1v46QLL8adVU9kyWk0fHgNTkVl2OGVGQFxcjvKmNVAjCkTqVSKp556il27dvG+d13Gso++ia4DB6h1otTbyrF5JYIt544lEGPKRvZ2s1//9i0AI0uxmzwTwYnZv09LocaUCVsIsXCEzETCXI5yZinUmBKiqnR2do650KEthFhAIrYnOpZAjCkZqsqGDRtobW2lvb2dpqYmWlpaWLFihW03W2g2CguwBGJMyejs7KS1tZVYLMZHr/wIAA8//DALFy60pBEC29LWEogxJSOVSvH6176Wt73hAoaeewqAi973bqIR+zMuOBHEtV0a7ZNnTImoqKjg7W9+E3u/9w/4A/0ARJ99iuZPXRdyZDORgNVAbBSWMaVizpw5SPsWJFhd13VdKlxh8MWNIUc2AwngurkdZcxqIMYUAU0Nouk0OtCLk6iFSAw56p+PiBBNVFNTU4OqIiI4jmOzzEMhEN5W5EXDEogxIdPUEKmXNzP4h1+Cl0bjCSrf+D763Arq6upGDcOteNVZRDc8ited2SbHnd1AxbJXhxX6zGZNWJZAjAmb+mkGN/wavDSpVIrk4T30P3YvB5vOZYeno4fpVlZx0gc/xdCuNgBiC5uReEWI0c9Mksfl3EuZvQLGhC2dgvQQvu+RTCbxPI90z0Eq4nFaW1vp7OwcuamI4FQkqFh6BhVLz8CpSJT9bOfiJBCJ5XaUMfvkGRO2SAxn9lwU8LzM5nIVTaex8fkXaG9vJ5Wy7WaLjgS1kByOcmZNWMaETGJxKt/wbmTjb6jc14Ezr5mhBady161/RVNTE9GozTcoPgKOdaJbAjEmZCIOUlVLxXkX4/Uc5tEn13PHP/8V1dXVthhi0bIEApZAjCmY4y2ECODE4lTPaeDc81Zy1tnn2GKIxSzPa2GJSBvQA3hAWlVXisgc4IdAM9AGfFBVD+btSfPAEogxBZDLQohgiyGWlPwPXvhjVe3Kun4N8LCq3igi1wTXv5TvJ50K60Q3pgCGF0Jsa2vD933a2tqOGWFlSog4SCSW0zEFlwO3BZdvA9491bDzzRKIMQWQSqVob2/HyZo7YCOsSpwjuR25UeABEXlaRNYEZXNVdQ9A8LPoqqbWhGVMnqiXBt8Dx0Xc0X9aNdXVrPv2N6murmHvnt38+w//k66uLhthVaoEJPdO9HoRWZ91fZ2qrjvqNheq6m4RaQQeFJEX8hLnNLMEYkwe6NAA6Y4X0N6DOLUNuAuXIdF45pzvU80QztbH6N3fQe3sefzFp/+U/UeSNsKqZE1oLawuVV15vBuo6u7g534RuQtYBewTkfmqukdE5gP7pxTyNLAmLGOmSFNDpLc+jb93O9p7EG/3i6Rf3oSmg+ap9BBDG3+NO9BDdXU1iVQvsfZnWNq0yEZYlTLHye04ARGpEpGa4cvA24BNwD3AVcHNrgLunqbfZNKsBmLMVAn4h/YB4PsequB3dRA5ZXiRQ0WThwFwgmYP7ekOI1KTJyKC5K/5cS5wV/BlIgL8h6reJyJPAXeKyCeAduAD+XrCfLEEYkw+xCpJ9R0ZWcsqPrsRHeinMhIFBKmahfYdGrm5U9cIVvsoXZK/iYSquh04Z4zyA8Bb8vIk08SasIyZKnFxTjmL/qEUnuch0TjSdCZ3//RnmWG6kRixc/44kzTcCE79QmKnv26kj8SUJlsLK8QaiIicDNwOzAN8MiMTvl4Ksy+NySauy2Cshti5b0OSPbiJGu574EF+eOd/8sY3vimz7Hf1bGKveWvmm6sqErMl2EubLWUC4TZhpYHPq+rvgw6kp0XkQeDjFPnsS2OOlhwc4qabbiKZTLJv3z6SySTNzc0jw3RFBCxplA9hOmail5zQXgFV3aOqvw8u9wCbgYWUwOxLY47W0NDAhz70IVSVgYEBmpubbSHEsiaZBJLLUcaKohNdRJqBc4EnOGr2ZTCxxpjQqPrgpfE8n/7BQQYGh6ivrz9mDasVK1awdu3acRdLNOVFrQkr/AQiItXAfwGfU9Ujuf7BBdP91wA0NTVNX4BmRlDfA1UQGTXDWD0PTR4itXc7/X09OLPmk4pUs3HjRs4++2xbCHEmsy8H4Y7CEpEomeTxfVX976B4XzDrkuPNvlTVdaq6UlVXWjOBmQr1Uuihvfj7d6CH9maWJBnmp/A6nqOvex+pZA+Du1+k0h9g07PP2kKIM5lYExaEmEAk89XtX4DNqnpT1qmin31pyod6afzu3WhvN6QG0N5u/O5dI0lE+4+gvo5sNQvgDhym/qQ5thDiDKaAipPTUc7CbMK6EPgo8KyIbAjKvgzcSJHPvjRlRAT6j4wuG+gZaZ6QeBUi4LruSBLxown6BwZtIcQZTcq+dpGL0BKIqv6WzGC4sRT17EtTRlQhEof04CtlbixTDhCJ4zYuJqE7SPb24NY14tc2sPRVno2wmumsEz38TnRjQiUOzpz5+F07R5Zid+YsGFkET9wIzF5AfPZ8Yp5PKp2mb2CQs846y0ZYzWQiZd88lQtLIGZGE8dBY1U485fjpwbBjZIcHCQRzxpdFeztIS7EY3HiiaqwwjXFxBKIJRBjEGHjs8/xk5/8hE2bNtHQ0DDmfuXGjGKfDVtM0ZjOzk5uvfVWnnjiCfr6+my/cpMDG8YLVgMxZmS/8my2X7k5EbUaiNVAjIlGo8esZtDU1GTDdM34RMCJ5HaUMUsgZsYb7vNobm7GcRxbCNHkpsSbsESkR0SOjHfk8hjlnR6NCagqnZ2dYy50aAshmkkp4uSQC1Ud3of9r4G9wL+RmZt3JVCTy2NYAjFlT1XZsGEDra2ttLe309TUdMwoK1sI0UyISDn1gVyiqq/Nun6LiDwBfPVEd8wphYrI3+VSZkwxUfVR38f30ixauIBTTjkF3/dtlJXJjxJvwsriiciVIuKKiCMiVwLeCe9F7n0gbx2j7O05h2fMNMskCw9vaDCzQKLvo+kU6cP7SXXvwek/widbPs6qVasAG2Vl8kAkt+OEDyOrRWSLiGwNdmEttI8AHwT2BccHgrITOm4Tloh8GvgzYImIbMw6VQP8blKhGpNnqppJFskeIPPVKVo9i3RvN/g+QmbJdj95mPe/9708+eSTNsrKTJGgeRhhJSIu8C0yX9I7gKdE5B5VfX7KD54jVW0jsxPshJ3oFfgP4F7gb8nsTT6sR1W7J/OExuSdKumBZFZB8K3P9wDBcR0SicwKug2NDTbKyuRJXvpAVgFbVXU7gIjcQeaf+bQnEBH5BpmV6cekqp890WMcN4Go6mHgMPDh4AkbgQqgWkSqVbX9ePc3plB8z0N9DxHBcRxUfcSNBvt6CNFohGhlFeloFWvXrrVRVmZKdGKd6PUisj7r+jpVXRdcXgjszDrXAWR3aE+n9Se+yfHlVAcTkcuAm4AFZHYIPAXYDJw51QCMmQpVpevAASpc6Dt8GNd1SSQSOOk0bs1JeH2H0HQKicZxq+qIOi6ViUTYYZtSpxznu/sxulR15TjnxspCuT/yFKjqbaMCEalS1b6JPEaunehfAS4AXlTVxWT267A+EBO6zs5Obr75Zg73DRCrriNakSAtEYhEETeCWz2HyKxG3KpZo/Y6N2ZqFD/H4wQ6gJOzri8Cdk9b2GMQkdeJyPNkKgWIyDki8u1c7ptrAkmp6gHAERFHVX8JrJhUtMbkUSqVYtOmTXz+C1/gB3f+iCf+8Aw/+vE9dHZ2AZnl2sVxEackhlOaEqGAp7kdJ/AUsExEFotIDLiCzLbehXQzcAlwAEBVnwHemMsdcx1GcEhEqoFHgO+LyH4gPfE4jcmv4XWs2trauO+++wBobm7msssuCzkyU+40Dw1NqpoWkauB+wEXaFXV56b+yBOOY+dRfYJ5nQdyOdAP/AVwH7ANsL9QEzpbx8qExdfcjhNR1Z+r6nJVXaqqN0x/5MfYKSKvB1REYiLyBYLmrBPJqQZyVMfKbePe0JgCs3WsTBgm1ode9D4FfJ3MiLAO4AHgM7nc8UQTCXsY+3USQFW1dmJxGjNxvp/5c/U8j4GBAaqrq0clCFvHyoQhH01YxUBVu8gsoDhhJ5oHktOKjMZMF9/3Sfb3c6SnFxFIVFSwe/duFixYYLUMEx7NDCEvZSLyRVX96ngTCqc8kdCYMKkqA0NDHOnpxUtnxmz0eUmqqqro7Oy0WocJzfAorBL3JTIr7m4DDk7mASyBmKLmez6e98qAP8/38YPDmDDl0kFe5PaJyCnAnwB/PJkHsARiiprrurhuZKQG4rqZOR2OzeswISv1JizgFjKjapcwelkTIVPJWnKiB7C/QlO0RIR4LEZtbQ1uJEI0FqUqUUmyr8+G6ZrQaY5HsVLVb6jq6WTmnizJOhar6gmTB1gNxBSB42036zhCorKSyngcz/MZGOi3DnQTOtWyGoX16cne1xKICVUu2806IuC6uK5LLGZ7eJjikMM6V2XPmrBMqDo7O2ltbaWtrc22mzUlI49rYZW0UBOIiLSKyH4R2ZRVNkdEHhSRl4Kfs8OM0UyvVCpFe/vobWVsu1lTCoabsU50lLOwayDfA1YfVXYN8LCqLgMeZvROiKbMDC+GmM22mzWlIE/LuZe0UBOIqj4CHL017uW8st7WbcC7CxmTKSxbDNGUrFIfhpUHxdiJPldV9wCo6p5gG91jiMgaYA1wzDdYUzpsMURTipSymEg4ZcWYQHIS7Cm8DmDlypX2Vhap4w3RHWaLIZqSo4pX7h0cOSjGBLJPROYHtY/5ZPZgNyVIVWlra+PFF19ky5Yt7Ny5k4985COjhugaU6osfxRnArkHuAq4Mfh5d7jhmMkaSqVI1M7itBXnc+aK83D8FLd+97ssXLjQahympGWasCyDhJpAROQHwEVAvYh0AH9JJnHcKSKfANqBD4QXoZksX5W+IZ+dB3pGyubUJLj0He+wIbqmLNh6niEnEFX98Din3lLQQEzeqSpHBlK4rovnZbZXPtiT5NRTT2OgPxlydMZMjWIz0SH8eSCmbAmu65BIJHBdF4BIJEI8FqO+vj7k2IyZIgUv6Eg/0TEVInK9iOwSkQ3BcWnWuWtFZKuIbBGRS6b8O01CMfaBmBLi+z5pz8dXxXEcXEdwHQcRmF0ZI+UpNTU1oMqsqjhuxLUOdFPyFCVduHVK/lFV/yG7QETOAK4AzgQWAA+JyHJV9QoVFFgCMVOgqvQPpXhpTze9/YPEohGWzjspkygch8pYhEWzHPpTHvGIS9TNJBhjykHI80AuB+5Q1UFgh4hsBVYBjxUyCGvCMpPWPzBIx4EeevsHARhKpdm29wBe8JflOkIs4lJXGaMi6o6ZPFQV31dSnk/a9/FtdpYpAVqgJqzA1SKyMVg7cHhtwIXAzqzbdARlBWU1EDNp4gh9A4OjyoZSaTzfB9ycHsNXpbtvcKQ5oCoeoToexbGaiilyE/iuUy8i2Tv+rQsmQgMgIg8B88a433Vkdg38GzL99n8DfA1oIbNr4NEK/u3LEoiZNC/tMauqkuTgK8NyqyriuDn2cagqfYPpUW3JfYNpKmMRnDH/PowpDgoTqS13qerKcR9L9eJcHkREvgv8NLjaAZycdXoRsDvXgPLFmrDMpFUlKpk/q4oFc2qpjEWpr61i+YI5RNzcPlaqjDR3ZRurzJhiU4jVeIPVOIa9Bxje+uIe4AoRiYvIYmAZ8OSUnmwSrAZiJk1EqIjHaKzxmD+7GlWfyngs51FWIlARdRlIeaPKojkmIGPCogqpwnzR+aqIrCBT6WkD/mfm+fU5EbkTeB5IA58p9AgssARijiPXhRCrqxKTenwRIR5xqUvE6B9M4zhCdUUUG+Vrip2iBRnwoaofPc65G4Abpj2I47AEYsY0vFf597//fWbNmkVDQwMXXHABZ599dl7ncTiOUBl1qYi4I9eNKQW2FpYlEDOOzs5OnnjiCb54zbUcGUgDMLsqzqFDh5k9e1Zen0tErNZhSsrwnugznSUQMybP83jvBz7Is3t66BkYAqCmIsaFyws+1NyY4qNWAwFLIGYcVVVV7OsdGkkeAP0pjwPJIRbG7GNjZjYF0jZa0BKIGVtNTQ0HhnqJuC5pzyPiulQmElZtNyZgNRBLIDPa8UZZiQjz6xJ0HKkZmZcRdSM0VsfDDNmYomB9IBmWQGYoVWXP3r08+rvf8bOf/Yzq6mpaWlpGbTcbdYXzFs1h95F+ABbUVhJ1rbfbGFRRa8KymegzUcrzOdA3yMtJYfFrXs/ffvXvqa2tpbW1lc7OzpHbuY5DIuayZE4VS+ZUkYi5uI59ZIwZroHkcpQz+28ww6gq3ckhfr2ti+d3dfGHnd08sfMwn/qzq9m1a9eY2806jtj8DGOO4qvmdJQza8KaYVKesu1AHwq4rks67dHZ0096Xg0XXHAB0Wg07BCNKXqqmZr8TGcJZAZyJbNrYCKRIJlM4nkesWiEyy67jIaGhrDDM6boZZqwyrt2kQtLIGVsrFFWsYjD8oZq9vcOQiRKdXUNC2ormF1bzYL62bbdrDE5UVs1GksgZWt4LavW1lba29tpamoaGWVVWxHlLcsb2dczQCLqMjsRI2Yr4BqTM9UJ7QdStuy/RonqT3kkh4Z3/ztWZ2cnd955J0uXLmX16tUMDQ2NjLJyHSERdVk8p4q5NRWWPIyZBM/XnI5yZjWQEpP2fbqTKX67/QDJlMeZ82o5vbGaiujoLWQdx+HPP/8F2g8NkFZ41/s+yC8euG/MUVbGmInJdKKXd3LIhSWQEjOUVn74h10kg02Ydh0eIOYKp8+twcnqv0jU1vGj3+9g/5EkAE9HXN57ydupzG2rcmPMcSi2cyZYE1bJ2d87OJI8hr2wv5fB9OimrK5kin51ibguQmZCU1uvT+2sWYUL1pgypZpb81W5JxmrgZSY2opj37K6igiRoyb6qUI0GsWtqSHzfUlwXQfBRlkZkw/lnhxyYTWQEpOIury2aTaOgO97nFQZ4fyTZx2TQObWxKmriOA4Do7jEo+6nN5Yg2szyo2ZsuEmLKuBFCkRWQ18HXCBW1X1xpBDKgoVUZcLTpnNOfOqOHikF9dPccft/8qqVatGLYQYdx0uPXUuOw4mGUz7vOqkKioj9n3BmLxQm0gIRZpARMQFvgW8FegAnhKRe1T1+XAjKw5HDh7gG9/4Bt3d3ezatQvP89i0aRNr166lsbERyCzHXhl1OaOxJuRojSk/ihZkKRMR+QBwPXA6sEpV12eduxb4BOABn1XV+4Py84DvAZXAz4E/V52ebFesX0lXAVtVdbuqDgF3AJeHHFPRSKVSbNy4kfb2djwv06He3t5uQ3SNKRDVgjVhbQLeCzySXSgiZwBXAGcCq4FvB1+8AW4B1gDLgmP1VIMYT7EmkIXAzqzrHUHZCBFZIyLrRWR99hLkM0E0GqWpqWlUWVNTky2EaEyBFKoPRFU3q+qWMU5dDtyhqoOqugPYCqwSkflArao+FtQ6bgfePaUgjqNYE8hYPb2j3glVXaeqK1V1ZbkuAJj2/DGXS2hoaKClpYXm5mYcx6G5uZmWlhZbCNGYAvI1twOoH/6yGxxr8vD0433JXhhcPrp8WhRlHwiZX/rkrOuLgN0hxVJwg2mPrr4hnmzrpqYiwvlNs6mJZ0ZUQaZ/Y8WKFaxdu3bM7WiNMdMr04SVcx9Il6quHO+kiDwEzBvj1HWqevd4dxsrrOOUT4tiTSBPActEZDGwi0xb30fCDalw9hwe4OZfvkRvso902uORLTX877ecRl2iYtSe5cMd5saYQsvfEF1VvXgSdxvvS3ZHcPno8mlRlE1YqpoGrgbuBzYDd6rqc+FGVRj9KY9fb+sKkkcaUNq7jrB970EOHToUdnjGGDI1kKG0n9MxTe4BrhCRePBFexnwpKruAXpE5ALJfNv8GDBeLWbKirUGgqr+nMwQtBnFEYg6Qjo9erkSBw0SijEmbIXaUEpE3gN8A2gAfiYiG1T1ElV9TkTuBJ4H0sBnVHX4n8aneWUY773BMS2KNoHMVPGIy0XL6nnspd309A8CcPr8WcypjDBNQ7mNMZNQiP1AVPUu4K5xzt0A3DBG+XrgrGkODbAEUpQaquJcf9k5/H5HJ5URZfGcKnoPHeDk5cvDDs0YwyuLKc50lkCKUDTiMqeqgteeMotUKoXneSxfvtxGWRlTJGw59wxLICHoHUyjqqgqqaEB/KFB6uvrRyUIEaG2tjbEKI0x41LwbEMpSyCF1jOQpvXxNn63dT9eapD3r1jE6TUpOjo6Ri2GaIwpXortiQ5FOoy3XHm+z2+2dfHo9i56+5L09g9y22PbiFbP5oEHHmCmLcliTMlSRloRTnSUM6uBFNBQ2mdrV29mFmuwCKICbQf68DzPFkM0pmQoajUQq4EUUjzicu6iWYiA62YWzow4wmnz60in07YYojElYrgJK5ejnFkNpIAcR3jNolm8/9yTefiFvTh+ivefM5+OHVu57LLLbDFEY0qFZmajz3SWQKaBqtLZ2TnmQodV8QjvOmseq09vJO156NAA6cZKWwzRmBLjF2BDqWJnCSTPVJUNGzbQ2tpKe3s7TU1NtLS0jN5uNprZoxyiUFURbsDGmIlTrA8E6wPJu87OTn7yk5/Q09OD7/u0tbXR2tpqI6yMKSMadKLncpQzq4HkWaSyhvd/7E+JR1327drJrd/5lm03a0wZKvPckBNLIHl0ODnEjfdt4fEXOkB9LntNE5/7P9dw6y3ftBFWxpSZcp/jkQtrwpqElOfTPzR6afW05/PA8/vYvKeHRCIB4vDj9S8zFKni6quvthFWxpSToA/EmrDMhBzsG+Knf9jF7oNJ3nzGPE5fWEt1RZSU57OzOwlANBqluqYGFI6kXc47bZmNsDKmjCg2CgssgUzIkf4Ua//zGV7YcwTf97n/mV18+fKz+OMz5lEZi/DG5Q388oX9ALiOQ0XU4ZyTZ1vyMKbcKKjlD0sgE9HdO8gLe46QTqVIJpN4nsePHtvKaXMrWNgwm9Pm1/K5ty7n3mf3kIi5XLHqFKri9hIbU37Kf5Z5Luy/2wRUxlxQn2QySTpYy6oi6tD+8svESNPY2MhFpzawsnk2glBbaR3nxpQjxTrRwTrRJ6QqHuHScxaOLIRYFY9w5etO4ZcP/HxkmK7rONRVxix5GFPOrBMdsBrIhFRXRGl502L+6FW17Ozq5dzmOfzu179g7969NkzXmBnGK0Anuoh8ALgeOB1YFex3jog0A5uBLcFNH1fVTwXnzgO+B1QCPwf+XKepumQJZILqqiqokyQbnv8NX2l9nEQiQUtLiw3TNWYG0cJ1om8C3gv88xjntqnqijHKbwHWAI+TSSCrgXunIzhLIEdRVQ4ePDiyGcxJJ510zFazy5YtY9asWVx66aXHLJZojJkJCrNZlKpuBnL+/yIi84FaVX0suH478G4sgUw/VWXfgcPsOzTAtt2HOG/ZXPp372XRgnnHJJHGxsYQIzXGhK0IOtEXi8gfgCPA/1XV3wALgY6s23QEZdPCEkiW7iN9/OC3O1h37/MAJOIRvv6nr6O2+gh1dXUhR2eMKRoTW423XkTWZ11fp6rrhq+IyEPAvDHud52q3j3OY+4BmlT1QNDn8WMRORMYq6oybZnOEkgWFZd/feiFkevJwTTrHtzC1/7kdSFGZYwpNsM7EuaoS1VXjvtYqhdP+PlVB4HB4PLTIrINWE6mxrEo66aLgN0TffxcWQLJMpT28XV0Au8b8BDH+jeMMVlU8dPhrbAtIg1At6p6IrIEWAZsV9VuEekRkQuAJ4CPAd+YrjhsHkiWqooobzq7CTfigoAbcXnfha+iNmGbPhljsimqXk7HVIjIe0SkA3gd8DMRuT849UZgo4g8A/wI+JSqdgfnPg3cCmwFtjFNHegQUg1kvLHNwblrgU8AHvBZVb1/zAeZpONtN1uXiPGVK8/jv5fWs2NfD286cx4XnNpILOrmMwRjTKlTUG9qySGnp1G9C7hrjPL/Av5rnPusB86a5tCA8JqwxhzbLCJnAFcAZwILgIdEZLlONY0HctludnZ1nI9dtJTBlE8iHsGx5itjzDEUfFtNMZQmLFXdrKpbxjh1OXCHqg6q6g4yVbBV+Xrezs5OWltbaWtrO+52s9GIS3Vl1JKHMWZMCqjv5XSUs2LrA1kI7My6Pu4YZhFZIyLrRWR9rvuNp1Ip2tvbR5XZdrPGmAnTwvSBFLtpa8Ka5NjmnMcwB+Oo1wGsXLkyp/F00WiUpqYm2traRsqamppsHStjzASFOwqrWExbApnM2GYyNY6Ts67ndQxzQ0MDn//853l+8xZe3LKZjo4OrrzySlvHyhgzIaqKX+bNU7kotnkg9wD/ISI3kelEXwY8ma8H7+lP8UInPNVVy6rXv4urTpvHrJpKW8fKGDNh5d6/kYuwhvG+h8zklgYyY5s3qOolqvqciNwJPA+kgc/kawRW30CKW+99ln9/aDMAdz+6jY++9XTWXHo2VRXWhGWMmQhFbU/b0EZh3aWqi1Q1rqpzVfWSrHM3qOpSVT1VVfM2AUZVueu3W0eV3fWbrfjhL4hmjCk1wTyQXI5yVmxNWNPGV4hHXfoGXun4isfcYlhR0xhTalRR60QvumG80yYWcfjTS1/NcHeHCHzy7a8mGrFZ5saYidECLWVS7GZMDaQiFuEdr13Ma5Y1smFbJyuWNjB/ThWVsRnzEhhj8khtJvrMSSAA1ZUxli2MsWzh7LBDMcaUMlUbhcUMSyDGGJMvVgOxBGKMMRNmEwkzLIEYY8xEhbyhVLGwBGKMMROmYDUQSyDGGDMZ1oluCcQYYybORmEBIOUwE1tEOoGXJ3i3eqBrGsLJp2KPsdjjg+KP0eKbuonGeIqqTmkJbhG5L3jeXHSp6uqpPF+xKosEMhkisl5VV4Ydx/EUe4zFHh8Uf4wW39SVQozlasYsZWKMMSa/LIEYY4yZlJmcQNaFHUAOij3GYo8Pij9Gi2/qSiHGsjRj+0CMMcZMzUyugRhjjJkCSyDGGGMmpewTiIh8QESeExFfRFZmlTeLSL+IbAiO72SdO09EnhWRrSLyTyLD21AVNsbg3LVBHFtE5JKs8oLGeFRM14vIrqzX7tITxVtoIrI6iGGriFwTVhzZRKQteM82iMj6oGyOiDwoIi8FPwu614CItIrIfhHZlFU2bkxhvL/jxFj0n8EZQVXL+gBOB04FfgWszCpvBjaNc58ngdcBAtwLvD2kGM8AngHiwGJgG+CGEeNR8V4PfGGM8nHjLfB77gbPvQSIBTGdUQSfxTag/qiyrwLXBJevAf6uwDG9EXhN9t/CeDGF9f6OE2NRfwZnylH2NRBV3ayqW3K9vYjMB2pV9THNfCJvB949XfHBcWO8HLhDVQdVdQewFVgVRow5GjPeEOJYBWxV1e2qOgTcEcRWjC4Hbgsu30aB30dVfQTozjGmUN7fcWIcT7F8BmeEsk8gJ7BYRP4gIr8WkT8KyhYCHVm36QjKwrAQ2DlGLMUQ49UisjFoXhhu4hgv3kIrljiOpsADIvK0iKwJyuaq6h6A4GdjaNG9YryYiu11LebP4IxQFospishDwLwxTl2nqnePc7c9QJOqHhCR84Afi8iZZJqEjjblsc6TjHG8WKYlxlFPfJx4gVuAvwme82+ArwEthYgrR8USx9EuVNXdItIIPCgiL4Qd0AQV0+ta7J/BGaEsEoiqXjyJ+wwCg8Hlp0VkG7CczDeWRVk3XQTsDiPGIJaTx4hlWmLMlmu8IvJd4KfB1fHiLbRiiWMUVd0d/NwvIneRaVrZJyLzVXVP0DS5P9QgM8aLqWheV1XdN3y5SD+DM8KMbcISkQYRcYPLS4BlwPagyt4jIhcEI5s+BoxXQ5hu9wBXiEhcRBYHMT4ZdozBP5Vh7wGGR8eMGW+h4sryFLBMRBaLSAy4IogtNCJSJSI1w5eBt5F53e4BrgpudhXhfdayjRdTsby/pfAZnBnC7sWf7oPMh6uDTG1jH3B/UP4+4DkyIzZ+D1yWdZ+VZD6Q24BvEszYL3SMwbnrgji2kDXSqtAxHhXvvwHPAhvJ/MHOP1G8IbzvlwIvBrFcVwSfwyXBZ+2Z4HN3XVB+EvAw8FLwc06B4/oBmebcVPAZ/MTxYgrj/R0nxqL/DM6Ew5YyMcYYMykztgnLGGPM1FgCMcYYMymWQIwxxkyKJRBjjDGTYgnEGGPMpFgCMSVPRHrDjsGYmcgSiDHGmEmxBGLKhmT8vYhsCvbd+FBQfpGI/EpEfiQiL4jI9wu5f4ox5aos1sIyJvBeYAVwDlAPPCUijwTnzgXOJLMu0u+AC4HfhhCjMWXDaiCmnLwB+IGqeppZbO/XwPnBuSdVtUNVfWADmQ3FjDFTYAnElJPjNUsNZl32sNq3MVNmCcSUk0eAD4mIKyINZLZCtZVYjZkm9i3MlJO7yOwT/wyZTYS+qKp7ReS0cMMypjzZarzGGGMmxZqwjDHGTIolEGOMMZNiCcQYY8ykWAIxxhgzKZZAjDHGTIolEGOMMZNiCcQYY8yk/H/MgjY71et5qAAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "ds_trajectory.plot.scatter(x='longitude', y='latitude', c='k', alpha=0.7);\n", "ds_selection.plot.scatter(x='lon', y='lat', hue='field', alpha=0.9);" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "You can also return the distances associated with each selected point by inputting the name of the distance variable you would like used as the variable named in the returned Dataset. Note that if a DataArray was originally input, a Dataset will still be returned in this case. Distances are returned in kilometers." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
<xarray.Dataset>\n",
+       "Dimensions:   (trajectory: 30)\n",
+       "Coordinates:\n",
+       "    lat       (trajectory) float64 -11.46 -8.57 -6.349 ... 37.13 38.05 41.56\n",
+       "    lon       (trajectory) float64 -150.6 -139.7 -129.9 ... 129.0 140.2 149.2\n",
+       "Dimensions without coordinates: trajectory\n",
+       "Data variables:\n",
+       "    field     (trajectory) float64 -162.1 -148.3 -136.3 ... 166.2 178.2 190.8\n",
+       "    distance  (trajectory) float64 176.9 33.55 71.99 198.1 ... 69.07 52.98 185.2
" + ], + "text/plain": [ + "\n", + "Dimensions: (trajectory: 30)\n", + "Coordinates:\n", + " lat (trajectory) float64 -11.46 -8.57 -6.349 ... 37.13 38.05 41.56\n", + " lon (trajectory) float64 -150.6 -139.7 -129.9 ... 129.0 140.2 149.2\n", + "Dimensions without coordinates: trajectory\n", + "Data variables:\n", + " field (trajectory) float64 -162.1 -148.3 -136.3 ... 166.2 178.2 190.8\n", + " distance (trajectory) float64 176.9 33.55 71.99 198.1 ... 69.07 52.98 185.2" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "ds_selection = ds_mesh.xoak.sel(\n", + " lat=ds_trajectory.latitude,\n", + " lon=ds_trajectory.longitude,\n", + " distances_name=\"distance\",\n", + ")\n", + "\n", + "ds_selection" + ] + }, { "cell_type": "markdown", "metadata": {}, @@ -144,9 +289,37 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 9, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "text/html": [ + "
<xarray.Dataset>\n",
+       "Dimensions:  (x: 10, y: 10)\n",
+       "Coordinates:\n",
+       "    lat      (x, y) float64 -26.81 52.21 54.89 81.46 ... -14.37 -6.231 56.52\n",
+       "    lon      (x, y) float64 -7.988 96.58 57.14 -122.8 ... 47.32 -177.5 -79.93\n",
+       "Dimensions without coordinates: x, y\n",
+       "Data variables:\n",
+       "    field    (x, y) float64 -34.79 148.8 112.0 -41.3 ... 32.96 -183.7 -23.41
" + ], + "text/plain": [ + "\n", + "Dimensions: (x: 10, y: 10)\n", + "Coordinates:\n", + " lat (x, y) float64 -26.81 52.21 54.89 81.46 ... -14.37 -6.231 56.52\n", + " lon (x, y) float64 -7.988 96.58 57.14 -122.8 ... 47.32 -177.5 -79.93\n", + "Dimensions without coordinates: x, y\n", + "Data variables:\n", + " field (x, y) float64 -34.79 148.8 112.0 -41.3 ... 32.96 -183.7 -23.41" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], "source": [ "ds_mesh2 = xr.Dataset({\n", " 'latitude': (('x', 'y'), np.random.uniform(-90, 90, size=(10, 10))),\n", @@ -163,27 +336,33 @@ }, { "cell_type": "code", - "execution_count": null, + "execution_count": 10, "metadata": {}, - "outputs": [], + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXsAAAEKCAYAAADzQPVvAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAAsTAAALEwEAmpwYAABef0lEQVR4nO2dd5xcZ3nvv89p03e2qvdiuci2bMvGBptmG2wgoSQhkEBowSGhJLm5IRBuyL03IZeQQOCShItpMQklJMbBCc2mYzDuVZZtyWqW1Xa1deppz/1jRtJKXmlnd2d2yp6vPu9HO+/MOec9U37nPc/7FFFVIiIiIiI6G6PZA4iIiIiIaDyR2EdEREQsACKxj4iIiFgARGIfERERsQCIxD4iIiJiARCJfURERMQCIBL7iIiIiCYiIr8vIo+KyDYR+YNGHScS+4iIiIgmISKbgbcDlwEXAq8QkY2NOFYk9hERERHN4xzgF6paUFUf+DHw6kYcyGrETueb/v5+XbNmTbOHERHRcdx3331Dqjowk23WxbNaCP2aXnvYK3xXVa+b1eCayEuvXqlHh0s1vfa+B4e2AZNffKOq3lj9+1HgQyLSBxSBlwH31nOsx+gIsV+zZg333tuQ9yciYkEjIntnuk0x9HnrwHk1vfb/HLinf8aDagGODpe46wevqum1Vu9nS6q6darnVHW7iPw1cDuQAx4CartSzpDIjBMREVF3DJGaWruiCqFqTW36fennVPViVX0+MAzsaMSYO2JmHxER0ToIYLavjteIElKfJJIiskhVj4jIKuA1wBV12fEpRGIfERFRVwTBbONZe61oncQeuLlqs/eAd6rqSL12PJlI7CMiIuqKCNhGZ4u9An6d0sOr6lV12dE0RGIfERFRdzrdjKNQNzPOfBGJfR0ZPjTC2OAE/Sv7SGQSiAgitPVCVETETBHAoPO/83U048wLkTdOnTiyb4hv3fg9epb2MFRw2b5viKeHc5S8oKYV+YiITsKosbUvlQXaWlqrEM3s64Cq8sAPHmb9xWuZ8ANGc5X4iWLJ45DA6p5UxZAZEbEAqHjjdPb3vWLGaS8isa8DQRAwenicDZdsoOCeiIcIw5BQwQsVs72nMRERM6LTbfbQfmIfSVAdsCyLzVedw+C+QRJO5fopgGVbGAvAMyEiYjLHvHFqae2KAoHW1lqFaGZfJ9aev4qnHtrNQCZOb3eKUBU3CEg4NtLht7RnIgjLhGEZMzBBDMS0EbP1v3bF8XEIAgLbwkkkcUyz2UNqGxZGUFX7zexb/1fXJiQzCc597tkU3IDhggsCPQkHQyo2/YVos3eDUXxvDHN4Anf3g4ia2KsuwFq6EbGdZg9vSjQMKR05wlM3fYGJvfvo2bKFFb/2WsJMmrhlN3t4bUO9Ju0ishL4IrCEir7eqKqfEJFe4F+BNcAe4LXHgpFE5P3A24AAeI+qfrc+ozmBAmELzdprIRL7OuKHyp6RPHnXBYVDY3nW96eZ0JBFmXTdj6eqBGGIFwbELLvlXDxz3m7Sbh+lR74Px7wSHv85RqoHs29ZU8d2Orxcjsc+/H8Y3vEkALm9e9DAZ/lvvoF4V7bJo2sP6hxB6wN/pKr3i0gGuE9EbgfeDHxfVT8sIu8D3gf8iYicC7wOOA9YBnxPRM5S1aBeA6oghG3mXhrZ7OvIeMnD84PjumaaBl4IMcsiVyzhevVJZqeqhGFAsVRkopDD9z1KvkfeLXM0l2coX2AkX2iqy2cQlvGDcYKjB2GS+5niEwzuQ8M6//bqROi6jO08OQ/V0XvvxfBbc7ytiEjFjFNLmw5VPaiq91f/ngC2A8uBVwI3VV92E/Cq6t+vBL6qqmVV3Q3spFIYpO6EVcGfrrUK0cy+jsQsk2Ofbcy2cEyTmGlS9AJMEfzAwzJNjDrc45ZKJSbcMopSDnzSCJ6CY5l4vodlWhTKLsmY05QZvyEOhpHA6Oo5qV8wMbIDIK05zzAdh1i2m8LoifQkyaXLFqQZbrYIYNf+fvWLyOT85JNzvZ+8X5E1wEXAXcBiVT0IlQuCiCyqvmw58ItJm+2v9tWVSrqE9vpOtOYvrk1JORZ9yTiGITimScI2idsmw3mXvUdzeKGSL5XnfBxVxQuCkyL4DCDt2JRyY/iFHMWJUdQr4weVu4kwDNEgT+jtJSz9HPWfQbVxS0wiQsZeh6ZS2GsuAsPCMOJYy87G7F/RsovWRizGpt/7PRKZDAIklyxl/dt+G6cr0+yhtRUzmNkPqerWSe10Qp8Gbgb+QFXHz3Doqb5Ydb/FrfjZGzW1ViGa2dcR0xCWZxMszsQB8IKQA6NFyl6AKuTKPkknPufjyBS5wB3bxiuX0PCEgBeLeZLJZKXPK+CXRhHTxoyfQ1j6Fkb8+WCtnPN4TkfM6seVcYy1W4itvRQQRAzEjjXsmHPFjMXouWQrl336M/i5CcxMF2YigdUGHkStQr2DqkTEpiL0X1LVr1e7D4vI0uqsfilwpNq/H5j8pV4BHKjbYCbRZuuzLXTZ6RAs0yBum6gqu4cmGC2U8IOKADuWiVWn6Kp4LI5tmCSdOAOZbkzTwokncJwTQmoaBqgSei5+/ihi2YiZICgHGPHrUbfx1b0cswvLyWI4SQwn0dJCfwwzHsfJZkkuX0GsqwvLjrxwZoohtbXpkMot4OeA7ar6sUlP3Qq8qfr3m4BvTOp/nYjERGQtsBG4u17nNWlkkc0+ooJjmfRn4hwZLxGGIem4zUAmjmPN3V9bRLAsi+50F6KK63n4QYBtmWQyWSYmxgh8j4TjYBgmfnECq2sp6hYJPRfDckCSYPRMf7CIiBlSWaCtm8g9D3gj8IiIPFjt+1Pgw8DXRORtwD7g1wBUdZuIfA14jIonzzvr74lzrFJV6wh5LURi3yBMQ1icSdCfihGEim0adZvVQ0XwRZWJQpGS7yEIKQTUI5POEHgudiyOhjnMWIKwOE5YLgIQuEVUFTN+Sd3GExExmXoFVanqHUxthwe4+jTbfAj4UH1GcHrqZY8XkT8EfpuKZegR4C2qWls18xkQmXEaiGkIjmWScKy6Cv0xFCj5fvVvJe+6lHwfMUzseBLCQSh+EbGcqtAfszIalceSrPuYIiIEwTKMmlq7ogi+GjW1MyEiy4H3AFtVdTNgUokTqDtNm9mLyCYqEXDHWAd8EOgG3g4MVvv/VFW/Nb+jaw9UFcsw8Ks+6xXvHMELAmK2jQY7wH8UtAyGBRpQmSQJGAaTJ0yqSm4kj2EKqWyqGaezYNHgKKgLEgdJIkbrr2ucEQFp47w3taL1s8dbQEJEPCBJgxaUmyb2qvoEsAVAREzgGeAW4C3A36nq3zZrbO1CGCrJRJxiqYQfBDi2TSzmnHBrNPoBRd2fYSYuJMjnQEwgjpXsOq71+bECj/7scWzH4uzL1uMVi5iOjRhmy7pIdgrq78OdGKM0LAw/9iQDF12C092Dk+5q9tBmjQDGAkiOU4/FV1V9RkT+lsq6QxG4TVVvm/OOp6BVbPZXA0+p6t5IXGrHsS0Krks8VvHt94KQIAyJx6ozQ3MZOFdC6ZtIsg+rez0EJmIlQCpukAA77t9FfiTHZdeeg3tgF+XQJ9nThdO/ArWcSPAbhIYFvPwRdn3jEe778MdBFbEsrvzbj7DimhdgtqsXkIB0uNhX/OznHjgmIj1Uon7XAqPAv4nIG1T1X+o4XKB1bPavA74y6fG7RORhEfl89c14FiJyg4jcKyL3Dg4OTvWSBUHScUjGHOK2TSYeIxU7YQIQI4skXoVk/hAxUogEGLE0YppI1V5aLpbZ9vMn2PrSC/FHDqOBDwpBqUQ4Pghhu+X2ayO0CKzkkX/4bMW9A1Df56FPfhY/V2ju2OaImFJTa19kJkFVZwocuwbYraqDquoBXwee24gRN13sRcQBfhn4t2rXp4D1VEw8B4GPTrWdqt547M0bGBiYj6G2JWJkEGs1Ym9CzN5nPW/aJj1LurEsE/W9SdsJ6run94OYJ8aPTjA2OE55fJygkD8paKztkTRg4RWKkzvxcsX2Ts8gghi1tXalkvVSamrTsA+4XESS1ZiCq6nk/6k7rWDGuR64X1UPAxz7H0BEPgP8V7MGthCwLIuLXryZYr6EnUgSFAuYtolpmUisud46g/uPMnpwiMXZkNyT92BZQvK8y7CWrGmL4KzpECOGkmPNy17O7ltvBQyQGOte/TIMp01NOFT97O3Oz/8fTONpUwuqepeI/DtwP5W4gAeAKVNGzJVWEPvXM8mEcywEuvrw1cCjTRnVAmLJmkVMDOdI9K1GJ4YgdDFjKYxMH2I050fruz7bf/EkWy5fSf7bX4AwwAUYPUj6JW/E7F3alHHVm1i2j60f+CP6LryAoQe3seyqy1nxoudhJxPNHtqcaOdZey0oQlAnw4iq/jnw53XZ2RloqtiLSBK4FvidSd0fEZEtVO6U9pzyXEQDEBG6+iqJvtRZWrEfizRN6AFKxTJiCOGh3TApHXLgB/j7HsfoXnx83aHdiXV3c9brX8P6V78cMxHHaPfzEkEWQNHldjMoNlXsVbUA9J3S98ZGHMstuZS8AMMySSdas0pSK9BMgZ9MMp0g05OCk1z+BdM0Mbr621bo3WCMcjCEVchQ3HOQsOTRdc5Z2Jk0dqozgtwWiuulNn/Jc0a0ghmn4YyNFTg0UebufcOYhsGVGxexJJvAsdrrw1pIGKbB6nNXYkpAfMP5lHdtw4nb2MvXYy3f0OzhzQovyHEw9236wst54E8+yNCdd2FKgsya9Vz2uY8TH+ibfiftgLRsuYK6cWyBtp3oeLEvTBQZK3p88c7dx2tG7jg8zruvOZteq/0X+ZqBhiFB6ONrSCCQsuaetnkqehZ3U8qXiT/npWQuuxpDQJw4Uoc00c3AC0dxwzHGH93J4J13ARCoS27vXp6++b9Y/9u/iWF1wE9SBKPDF2hVZdpUCK1GB3yzzkwYhGw7PHFSceBSOeDxwzmeu25+xL5UKFMYK2BYJrGkQyLVnmIFoIGPP3GUcn4YLAcju4hB32Ug3piIz3gqBnTGRTnExxALd2R0Um/li+kOHUWDsCN+kbIAgqqgfonQ5ov2Gu0sEEPIxqxn9XUn5udXNTGSY/ej+zDjDgU/wAuUUsmdl2PXG9WQIHeUwvDTuOUcbn4Yd3APMTEYc9s7CGg+cIwsoboMXHEZTnc3AIbYiGWx4jUvx4x1zlqSYRg1tXZGtbbWKnTAPOLMJNJxNq8wePTgOLsHc4jAuat6Wds3P8m+hg+N0r92MY8eGMULQmQwx7rlPSx1rJbP+heEAQW/TCEoEyIM2Em0ME4wqUZP6JUwAh/f7Ozb9nrgmD2syLyKYvEAz/3ip9n7xa8TFl1Wv+41pFY3rmLYvCPtHh07PcfKErYTHS/2hmGQjFn8xuVrGCt6mKZBV8Im4czPqcdScQ5OlPGq1aoUeHpogiXdre95EYYBlAukfI/AiRFoiGnHMIpCeEzwxcAybUJpoSlMC5OwlhJLDxCmPM794/egoWKnOyvLqABGh/vZ02JVqGqh48UewLIt0jZNcblMdSdxD51cH9kwTQJVWjlGMgh83LEjuKWJSkcBguxirO4lJMp5iuU8oUCsZxkeIVk73dwBtxGGWBhiVZLZdiILxmbfXue4IMS+mViWyeLeFPmDY5XHjkU2HcNq8dwnguKVcyf1uYUxzFgKZ8lGrMAFw8IjJG5GmTEjJiFg2O1l4pgpSn3SJcwnkdg3mEQyxhLHwnFsRoseCcdkSVe8IZWr6okgWGISqH+8zxCplEM0LUyz8tWJLPURz2aBRNBGfvYRp2JbJgMZg75UZQZstoM9Uwxi6V7IDeNrgCAkMv0YLRJhG9G6iHS+zV6rKY7biUjs5wlDpCEh5KoeaIAYs/fd16p/2GRTjBgGVqobM5Yi9F0MJ1GpXNXiHkQRrUGnJ0KDupYlnBcisW9TwjAAlCAoYxgG4g8hRhIxal/1U9VKNKxX8fs3bKeyr6roi2EijonRphGrEU1CBCMy47Qckdi3KUGoDOeGCas29WwiSdzKz1DsQ8r5Ccqejx+EGIZBprsbEGwrMtdEzA6h871xVCOxj5gHwlCZKBWPCz3AeLFIrCuLhoWaBF9VCV2Xsufj+ZUUwkEQUCgUwbTIWu2dT30yGoZouUh575OI7eAsX4sRb2+/x/zoGLmDhzn65C5WXLoFpytDvFX89QWkwycL9cxnP19EYt+GKJWZ/cl9lX8z8Y+p7Cd8Vt9Evkw21TliH+THGfrXf6Q8OooIJJevovdVb21bwc+NjfPQF77KnR//DABmzOFVn/koq6+8DKMVIpkXihmnTjZ7EdkE/OukrnXAB1X143U5QJXO/0Q6EBGI2Q6Thd0ybUSsmhdqRQTTcbAmiYNlWYhp0UlrsBoE5B++i/zRo/h+gOcFTOzbQ3n/rmYPbfb4Pvf8vy8efxiUXe76+I0UR8fPsNH8IoZRU2tn6lSDFlV9QlW3qOoW4BKgANxS7/FGM/s2xBAhHYshdFP2XSxDSMfjCDNLWSBikOrqRopFAlUM02LH/iHOXbOk5n2oKsMjeUplj2w6hhOzcVqofqrr+vjlEpPfmjBUgnJ5zvv2cjncoUHGHn6Q5LoNJFatJtbVmOyfkwmDAL98cjK98kS+4cetFVkAlaqUSprjBnA18JSq7q33jjv7E+lgDENIxR16Umm6kiks08IwZiayIoJpGohp4YZQKntsXruEZKz2/TxzcJT9zwzTFRckd4Rw9BChW0LDFinaZhjEzrkUsU+kyjAzWZzVZ81pt6HnMXbf3Tzy3v/G3n++icHvfhPcMmMP3MvYow/h53LT72SWGJbFhquvOqlv82t/mdgUla4mSmWG8wX8ef48xJCaWvsiqNbWgH4RuXdSu+EMO34dk2py15Nm16DdA0wAAeCr6lYR6aViv1pDpQbta1V1pFljbGUMEeZqNhQRMsk4meTM3Std1+Oxxw9wxcUr8Q/vOp7PVco5nKUbaAV7UCxmk7eT9P36O3EfvxexY8TPuxQPi7lkSvKLRQ7eeguqSmLxYla89vU8+YE/JsjlKIVFshvPY+N7/wwrnanbuRwj3dvDtR/5M5bf/E2OPvYEa699AauuvAwrXsn7HwYB5Ykc5UKRwLJ4ZGiYcT/gqg1r6JmPQuYiGHbr3N01AgX82n98Q6q6dboXiYgD/DLw/jkM7bS0ghnnRao6NOnx+4Dvq+qHReR91cd/0pyhRZyJcjmgqyuO6Y7jT0rcrUGAliaQdG8TR3eCTHeGkVGT4obnYjkmtp0gm56b6AkQehVvqL7nPIeRn91B6Zn92F3dqCpjO7Yxvu1hep/zvDqcwbNJ9fVw0Vt+Ha9UJp5OnRQQN7bvAP/5e+9n70PbSfb38ZK/fj/ji/v5+e69XHfOWZgNvghXsl42/0LfaLT+uXGuB+5X1cP13jG0phnnlcBN1b9vAl7VvKFEnIlk0qEnmySctFAsUqkfi9EK84gKtmWyqD/D6jWLWL6sj2x27l44ZiLOkutfhgBmPI6fm8BwHIJqtv9QlWBiYs7HOROWbZPIpE8S+tLoGD/5y09w+LGdqEJ+8Ci3vfevuHzZYh4/NEiuDmsV0yKAKbW1Niassc2A19MgEw40f2avwG0iosCnVfVGYLGqHgRQ1YMismiqDat2rxsAVq1aNV/jjZiEaRosXdKNbQlaHEO9Mo5tYsSTSLxFfL4bhGE79F71IuLLV1E+sJ/uy65g+I4fUypXPGKcTJbsJZfO+7g0VIZ37KmY+Krkh4bxiyUG0mkccx5+8gthgbbOQVUikgSuBX6nbjs9hWaL/fNU9UBV0G8Xkcdr3bB6YbgRYOvWrVHljCbRlUkQ+AHW8vXgVVI3YMeQDk+YViq55EpKauVq0hs3EgYB6//8f3Pktm9h2A5Lrv9lzAbY66fDsG1WPu8SRp8+gG2ZeEHAwMZ1EHO4cv1qEvPkKSXS2WJfSQJev3NU1QLQV7cdTkFTxV5VD1T/PyIitwCXAYdFZGl1Vr8UONLMMUZMj2mZgAlWZy/KHcPzA57aPciyhM+ef76R/K4d9GzcyMob3sOKN74FNQzi8eYUc4llUjz3j98BwL477qF7/Rqu+h/vwe7KsHQKb53G0Pkze4CwzaaYTRN7EUkBhqpOVP9+CfC/gVuBNwEfrv7/jWaNMaL5qO+ibgktTGB09YJpI/NhijgD4+NFuhzl8Fc/T37XDgBGd+7E+OwnWfXu92E2SeiPkezt4fkfeA+B64FhkOhuvO//ZETAtDv7zq6yLtNeF7Rm/moWA7dUF5cs4Muq+h0RuQf4moi8DdgH/FoTxxjRRNT38A/sxNv7GKqKYVk45z4PIzvQ1MpYYahks0kG9+4+MVZVyvv3toS7KYDTzDw5Im0fHVsLUVnCGlHVXcCFU/QfpRJF1nK4fkje8zkwVmJxOkYmbhHr8IRPzUTDgMKux/DKZRQwxIM9jxI/73lgx5o2rnjcZmJogvS6jYw/9ghQKT8ZX7MeWiWYrMl0vhlHGhVB2zCavUDbNvhByFNH83x72wE8r1K56epzlrB5SYaYE72NjUCD8LjQQ8Wd0S0WcUJtajnETDpOEHTT9Ybfxvjy5yjs3kl23XqW/ubbMBKdk0Bu1kjni70yY7fKphOpVI14ofLzp4Yol46lFVZ++uRhzl1SP4+LQq5EseghpoFpGsRsk3hyLnGe7Y0XKFbfEryjh473mYvWUPaVZPMm9gB0Z5N4qRgr3/4eLMsEVYxkKiq8DlTUvj5iLyKfB14BHFHVzdW+00bZi8j7gbdRicp/j6p+ty4DmYJ2m9l39uW3nqhSdL2Tulw/JKzTkrxb9vBDZdgLeGIox67hPOUwxPP86TfuUEqhYG64lNja84ktWU3s7OcQ9K7EjrXGBdC2TGJdXZjJFGYqHQl9lWOBdbW0Gvgn4LpT+o5F2W8Evl99jIicSyW3zHnVbf5RRBpzE6jMJDdOSxCJfY0YwPnLuk/q27wsS7noTvn6meJ5AUfyZXYP5pgoegyOl3j46dG2q3N5JkLfJ3BL+IVxQt+bNllaNpNgrBDwVKGLp62V7M3FUMNqySpaGgT4+QmCUrHZQ2k+IhiOVVObDlX9CTB8SvfpouxfCXxVVcuquhvYScWdu+4oEKjU1FqFyIxTIzHH4pKV3WQdk/1jJRanHNb1JXHqZJsU02A4f/KFo+j6lP0ApwPc2HzPIxg7QmHoACjYjkNixVmYsdPbuEWEFct66OtJUXZ90qk4ltV685OgkCf34N1MPHA3Vk8ffS99JVZP34LwSJkKYUZBVf0icu+kxzdWAybPxOmi7JcDv5j0uv3VvoYQeeN0MKm4zYbeJMtTDgZKzDRId9VnQc40hHTCYXSS4NuOidOCs9hZoSH5wYMcSyzvumWs4QPIotUY0/jNJxIOiURrmG5OJQwCJu7/BYf/vTLRFITS7h2s/IM/w0w119++aYjMpAZtTRkhaz3yFH0NCn1qLRNNLURiPwNEhHRXgkb8hGMxm7UDafKuz0TBw7YN1i/qwmzzZFHHCHyfU/0XAs8jDEPaObOCX8gxfO/PyPsVryFbTHTwIN7I0QUs9g33xjldlP1+YOWk160ADjRqEG0WQBuJfSsRt022rOrFC0Js00AAqwNMAaoBhmVjxhIE5RM2bTPT2/Y5dIoEaFcXYTXFs6s+hmk1JS9OK9FgE9bpouxvBb4sIh8DlgEbgbsbMYBKBG17TcQisW8hRATbFOwO8lHW4Cjq/hw/XE982Xq80SOo52KmezDTWcwGFsjWwIfAA8NCGpS358nyMOuvexX53U/iDg+BYdBz7SvAaU2z07wggpj1eb9F5CvAC6nY9vcDf05F5J8VZa+q20Tka8BjgA+8U1WDugxkCiKxj1iQ5EbziFRmdLG4g2mbaFhAi98EsYlZASpD0N2D4oBhYTewmpG6Jfy92wiGnsHI9mGvvwiJzT0RWFAsUNr/FP7QEZKbLuCsRD/fHX+Ua//of+EOHcbp6sZ1bIzEfCUda0WkbmkjVPX1p3lqyih7Vf0Q8KG6HHwa2s1TLhL7iDkzNjiGHbMIywUQ8I00YGMYOSR2FerG0dwEGD5OykdNxTAbF2mqvov35L0UdjxI4IfIgb0kxkeIXXwNhjPz8ovHCIp5hr75Fca23YWimD+4mWW/+Qec07OKj+/9IX2xDEsK47y463zMjk/xewYEaOAdWytQ73z280Ek9i2KH4SECpYpJxWimA2qSqhQ9nxUIWabmIbUJQioXHRxEg6l4SOEXjXobHyM7PIVqJEGTwnH9h9/fVDOYQ6cM+fjngkNAopP78AtVwPSAsjv24WzZW5LakExz+i2X6DVpTnfLzPys2+z6VVv5p2bXoqqkrRiZOyFnjJBoEGxTK2DQH2Ll3QDnwU2U1kSeKuq3lm3AxCJfcvhhy6uD7uHc5SDkMXpJIvS8TnZ8UNVDo3kyeddQlUSMYtl/Zm6BCd5ZQ8N3BNCD6AQlAqYmS7C0ill8NRAfbduNt2pCIMQtZPA2InD2jFU5yb2GgbHhf7EsTwUZVG8e077rjfFfAmv5BH4AalsEic+z2sIHeBYMB11zmf/CeA7qvqr1cLjdbcDdv4n0mYEofLAM4d4ZnyEofwY248McjhXmlNahkLZJ5crE4Qhqkqh5DGeL+OWvek3noZY0qkWL5nEcdc7A8w4GGmQBBgpkFhDhR7Aw8I553LEqSbQMU3i5z2XUnluqavMTDdLXvZG4svXVjrEIHPZiyA2e9NQIygOjzL66OM8890fEuYmePqxfeTH8vN2fKmmOK6ltTOK1NSmQ0S6gOcDnwNQVVdVR+s93mhm30J4QZ6C55OfNEsO1efwRJFF6TjOLBeEwlCPuwZO7vO8ACc2N+G1HRsRCNIp3HweESGWjGPF04hhYKR7CNxixSsGkGSWRjvWJ5Ixjo53kXzh6/HHRzDTWcbGy/RUF4RDt3TinTTMmjx1wmKe8s5HMctFlv3Sm/GLOXxLsHoXE7NaR+zLYxNs+8Rn2faZSt1qK5Xk+Z/7GOPDXaSy85XjXsDqbG8kBcKw5t/jdFHC64BB4AsiciFwH/D7qlrXK3Qk9i2E4uOYJsLJARtxy5jTLVjCsTANg6Cai8YQIWGbGEZ9bI6WbZPqX0Qi64KGmLH4iayHhoXZt6Ii9mKCYcyLb313fxejR/NMuEniJaVnUQ+JpINfylN84PuUD+4kFksRP+9KrOUbEfv04uQX8wzd8nny+3ZjGoJz/0/o+5UbMJcswTFaR+gBgkKBJ//55uOP/XyBJz79RTZ/8I/nbxDCgkgKN4N77emihC3gYuDdqnqXiHyCSnK3P5vTAKc4SESLYBtdeMFR1vZ2sXt4HAUSlsPa3jTWHGz2AqxclGEs74IqyZhF4AWke+o30zNMc0oPG5HqYt08B0+Zlknf4i76Fp8oyTdRKhDsfIDi05W69l4wRvjg98ksWnVasQ9DpXR0kLHdTwHgB+AFJeIP/Izua34FibeWqAXlMqdew91cnkRqPnNCy7x/3s2gjq6X+4H9qnpX9fG/U83kWU8isW8hRATHirMs67K0azFuEJKyE1hzrLlq2yYi0BW3UVU01LoKfbtQKhcxR07kxkfB9T3CiWGM5NQRr/lSmSB8tglMoZLLt04EpRLqeUAlL74xS9fFWE83S664mAN33IvvBSCw6Q2vIdaTxQ3GgQDbyM4kUdksWABiX8f0xap6SESeFpFNqvoElRiCx+qy80lEYt9i2EYa24AgLJGwLQypz0dkWWalyMYCpgB0DayEg7uO94llYWT7T7uNiKCpLInVGyju3VndxiZ58VUYdVqY9XM5Dn3j6wz/7A6sTIblv/4bpM89DzM+8/072S4u/79/ye5/+09yu/ax/PoX03fRZvLmkwyPVyaOaWcdffHLsM0GpXRYIJWq5ujcdSrvBr5U9cTZBbylrnuniWIvIiuBLwJLqGTIulFVPyEi/xN4O5UFC4A/VdVvNWeUzcNsMVtwJ5C04+jys0gUc5SffhwzniJz/vORMywmphMxdo3k6Ln+N4g/vRNyoyTPuhCpU+6bMAg4+rM72HPzzQRhiMghch/9G7b8/T/OSuwB4r09bHrrbxC6LlYqScF7Gh3zWORejVcsYHfHKVj7yJrn1eUcpmQBBJXVU+tV9UGgXtk/p6SZM3sf+CNVvV9EMsB9InJ79bm/U9W/beLYIjqQgWSKwSKU1m6he9NlWGIgTgyZxky2clE3g2M5xnpX07ViE2EyRjxZn4uxVygw/MB9+EFl8VxVKebz5J56it6tvbPer2FbGHblvDRvsPsz3+OJL30N9X36LzyfKz/5YYL+MqbRAFu+GGe8gHYKYdheF7SmjVZVD6rq/dW/J4DtNLDQQMTpOSY0C4GBRIpl3X0kkxmcRAq7hvUQ2zJZ1pflnFWLWd6fpatOQg/gGybxtetP7hSD2IqVU28wC8qHJ3j8pi+jfiWieOihR3jqK7eiXgOT9BpSW2tntMbWIrTEpUlE1gAXAcdWo98lIg+LyOdFpOc029wgIveKyL2Dg4NTvSRiGlw/wPUDyp7PnkPDHB2fv8CbiAqmZdP74mvpv+K5iGFgZTKsedvbIV6/lAuFpw8jk9Z+BJPc3kOEboPqGwuIYdbU2pl6BVXNF01foBWRNHAz8AeqOi4inwL+gso18S+AjwJvPXW7alDCjQBbt26d9fXTDwIOjeTYvvcwyZjNppWL6J+34JPm4YchQxNFRvJlTEPoTcUYGstjW2ZdZ64RZyYesxm3HDKv+y1Wv/13CIGJso+ZrF+0fP9F55PIDODmxqs9BquvvwY73ajveefnxlGt+wJtw2mq2IuITUXov6SqXwdQ1cOTnv8M8F+NHMPewyN84+fbuPSc1Szt7+bIeAHHMulKnSx4BddnpOiRK3usyCaJ2wZmm4Z7h6pMFD3GCpUSiEGoDE2UWN7bRb7kRmI/zyzqyTBqWxQ8H9MwWJzNnjFvkZ8vUB4cZOTue0hv2EB64waszOmrYlmpFC/43N+x7R//CXdsgnW/8jIWXX5xYwOf2vS30ck00xtHqOSC2K6qH5vUv/RYMWHg1cCjjRqD6/k8vOsgv/aii7n5wQN85u7HWd6T5IYrE6QS4XExz5U9vv7wM/xo12FUYUk6wR+9aBP98xqoUj9UFdc/2U6vQNkLOqpwSjvRna7NbBMGASP33Mvj/+cjEIQgsPyVv8yqN78BKz214FvxGH3nn8NzPvwBNAixMnHMBhVzgWpunAbWKmgV2q0GbTN/2c8D3gi8WEQerLaXAR8RkUdE5GHgRcAfNmoAoSrrlg/wHw8f5NvbDnF4osz9+0b4yO1PUHBPFLgZLXl8/6mDBBoSEnIgl+c72w9R9htWBKehiAjxU4qZGyKk4zYJp/N/pO2MP5Fj71f+lXy5xLhXJu+57P+vbxL6NdjfUx7Dzs/ZNfE5DuS/Qzk42phBSjWoqpbWtlSCqmpprULTZvaqegdTV4OfN5/6uGOztC/LQ3fsO94nAkdyLgUvIBOvCN9wsfysbUeKZbwgINaGgUqGCOmYTZAOmShWat32pmOICMn5ToXbZvhhiBv4OKbVlPrAgYaUXJdAw+pjpeiV0WmyovphidLoIE4wgB0b5lDuOxT9g6zOvB7brL/tfkHkxmkzm/2Cv2eP2yYbFnVhmwaOZZKMOWTjNkn7hIiv6k7SlzxZBK9Y20fCbvr69qwxDaE7EWNFb5rlvWlSMZvEHDNgdjrFYp7y2BClHQ9RPnqIcnH+vZdcx2bpq375pL7FL34RgXl6cfXzBYbvfJBH3vdxtr//MzgPLmaV/XrGSo8SaK4Bo1wIM/v2o33Vqk50JRzecsUajhZc9hwt0J2wefvz1uJYJ66Dlqn8txdu4rbHDzNe8rh8TS8b+9NzWqBVVQiPgpZAbDB6qaxXzx/1ynq5ECi6ZdyndzD8o3+HavbQnsuvh3O2EovNX2UqTwS97CIu+NjfMP7zXxDfuIGurRdTjsdOW+2ieHiQu9/9fjx3DFCO3nc/V331M3QtPoupb67niLAgImjrWalqPljwYg/Qn47xwevPxQtCHNPAtuQk80zajgMlXnn+EgJVkpZJ2pnb4qz6B3FDE5UMQogVPINpLUUaEdHYRLRaMCX03UoeGsNscBKuxmAFHkfv/+FxoQcYu/+HZM6+ZF7H0RNP8ExhnJ8nhcx1zwcxOM8QzjrDBefQ7T8GP8QQh1DLaBBw8Ds/YeXvXoNtZBswSul8sY9cL9sT0xAy8TO/FWk7TrpOE28Nc+SDND946jCDuUMkHZsr1wywNONid5DYqyqB71IqjOGGAQIkk1lMOzbnTJ7zjSCo557Up75f18yXtWAaBuf3LWFZqosgDLENk95pArDS69cAgiExDCxCArIbN5Fy1jcmXQKgHW6iqSRCa6+ZfYdfflsTNzC5++kRBnMFAAqux092HSHQzhF6AFRxS3lKgUeoIYGG5Aqj7bl4Z5p0nXNynqrspi1zrms7GwwRBhIplqQy9MTi5CdKlIruaV/fd+lFLHnxlQiCiMXSF7yQxVc9F6uRyfZEamvtTJulS2iv6VWHEGIxXDi5/mvJD3EDiHfQGqkCXnCyS6CqEmiISXvN/CwnTmrzFdjdA7j7d2IvXkV87bmY8brXha6Z0ugYYaAUxwvs3jvK0lX9LF/Tj3HKWpKTzXDB/3ov5/z334NQsbMZnGzXafZaB2QBmHGgpYS8FiKxbwKmYbKkK8tYyaWS/NMgE0/jzLJgRSvjxBJ4xRMXNsdJkvfLOA0uOt4IrHgSc/1mEqvPBtNqasFsb3SU3f/w/zj00ztJLl/G6hvezo7tz5BMx+lb9Gwhd7JdjRX4SSigC0Hs24xI7OtE2Qso+iGHxoss7YoTt83T+uA7psHWFb0EITwzXqQ7bnP56j6cDoteNQwDK5YkLqCBT2gYYNkYLZQcaqaIGHCGerXzQVAssv8r/8rB239IGIRMPLmTJ//qw5z/D5/k8KGxKcV+flkgM/s6fo9FZA8wAQSAP03N2lkRiX0d8IKQbYfG+dI9+yiVfGxTePMVa9m0KE3iNHaZhG1yxeq+4zbfmGW0py17GmzTomzZeIYQqOIgZJ0Tpg9VZaRQxA0CTMOgJxHH6sA7nHoSui4Tjz12ksnbHR7BGxkhc4aqW/NKhy/QolRKLtWXF6nqUN33WmUhXH4bjhcqtz5ygHzeJQhCSm7A1+9/+ln5Z04lZhnEbZO4bXak0B8jbSfojWUYiHedJPQAR/N5TIGumEPJddkxOIwXtGcaivnCsG1SG9bjxOzj3xsrkyG5aIBsd/PWEI4jgopRU2tvpMbWGkQz+zpgAKP5k70hxkt++znizjMl1yUhSml8CDfwycQSdGe7GS+V6Uu1gGi1KGYyyco3vgH36DAj9z+A3dvHunf+LmYiTqKrRd63thfy6ZnBz7tfRO6d9PjGaor2k3YH3CYiCnx6iufnTCT2dUBVuWhlD3fuOFFE5ZJVPUyMFclmandvU1X8UDENwejgmf4xbNMkf3SYMKzM5L1yEcMwiSWbbXNufZzeXjb+9z8C04BQMeIxjFbKNLkAvr8z8MYZqsEG/zxVPSAii4DbReRxVf3JnMZ3CpHY14GEY/GaC5fTl7DZM5RnbV+SS1d048zAbFn2A/aNFnlmrMTSTIzVvUnibZhkbUZoiIGeZPoMfJdEZLOvCaurPkXP689CWKCVuqZLUNUD1f+PiMgtwGVAJPatSCZh88IN/RRXdiN+gGMZdNVoP3WDkAcOjPHggVEUn21HfM7P93LpiiwJe/7yrsw3hmESj8UolssEYYghBslE6ll+4guJUENCDbGM9v5p6gKY2dfrDEUkBRiqOlH9+yXA/67T7o/T3t+oFiOdipGeRUETVeXxwxOEWibUiu3/scODXLoiixeMY5sdatYQiGX6MIxRNPAx7Dh2qqvta5POBlXlqHuEx8cfohgU2Zg+l6WJlcTMOGEY4I7nEMOglC8hiVTNE4mmIAJtfrGalvp64ywGbqkutlvAl1X1O3Xbe5UO/0TaAwUsE0LvRPCRbRoE6qJhvmPFXqTid+90nXAXbGagUjMZ80a47dB/kPMrdWL35HfwokUvY31yE6X9+9j9mf9Hbs8eus49l2VvfDtjYUi29/SlCE8lDD2KwUFGyw9jSZJsbDNxa1GjTqftzTgiMsEZrPJ9G9bX7Viqugu4sG47PA3t/Yl0CJZhcMny7pPcLy9Z3st4+YGOdsmEagk7wzjeFipHygfI+ePEjAT9scXEzARPTGzDLxbY86lPMLHjSdRzGXvoQQ599Sb8/Mxy6ee8Xewc+RSD+Z/y9MTN7By9kZI/OP2Gs0WM2lqLoqoZVe0CPg68D1gOrAD+BPjLyotqbC1CNLNvASxDWNeXYlFmLQfGx1mcjmEYh8l5h7Hl0mYPL2IeMMXiOZnnsj62AXf4ME7vYoZ1BIKA4jPPnPTa3I4nWT5NltbJBGGZoFjmLOfdqOfjGRMc5Dvk3J3ErYF6n0rVz75jJikvVdXnTHr8KRG5C2gpIa+Fmr4xIvLXqvon0/VFzJ6YZWEaBglbmXAfxTC6WJa+DqsBJeMiWo9lzgrKT+fZ96X3E5aKGLEEK3/jnchqk9SKFYzt3ns82jq9YSNiGHh+gF2Dx1aYLzP+nUd46AtfICgW6L30Ms760z+kKIcbd0ItPGufIYGI/CbwVSry/noqKQ1opYCpWqj1E7l2ir7r6zmQUxGR60TkCRHZKSLva+SxWgXLSJK0l7M49QJ64xfhmI0oLBHRijg+jH/rFiw3xDYcHB9G/uvfENNk1Q3vovusDViORe+FF7LmLb+NaEB5+Cil8fFp9+1P5Nnz6ZsIipWU2sP33M2hW75DRjY07oTqlOK4BXTgN4DXAoer7deqfZ1lxhGR3wV+D1gnIg9PeioD/KxRgxIRE/gHKheZ/cA9InKrqj7WqGNGRDQVMQjGRrGME4FR/lgl939i1RrWv/eDYBiEnsvB7/wXz/zHzWAYLH3xtax83RuwUqdfrC0+vR9RA0tSFW8vEYq7nsHwbGhISntB6+CN0wo6oKp7gFee2t+/cUNLCXktTPeJfBn4NvB/qCxSHGNCVYcbNqpKQMHO6io1IvJVKm94JPYRHUvq3C3kHrzr+OPkORegYYhhmljpNIeOjlG+8ycVoQcIQw7/4Db6r3wBmU3nnHa/6Y0bMBNJgmIRUyo/+b4rLsdMNtJ9sy4mjqbpgIh8kjPIef/69cg8i72IvAv4kqqOzGb7M4q9qo4BY1TsVFRDeeNAWkTSqrpvNgetgeXA05Me7wcmL5IgIjcANwCsWrWqQcOIiJgfzESSvl/6dcxsL6U9O4ivWkv3C67HTJ5Ys4lbBvlDJy/WigjuwWfgDGJvJhOc91d/zp7P3oQ3Ns6ia17IwAufj9GgSGWd2QLtmfLGTKsDDeTe6V8y7zb7JVTubu4HPg98V2dQKq3WBdpfAj4GLAOOAKuB7cB5Mx5ubUz1Lp50UtUvxI0AW7dubdg1NghDcuXKekzKMbE6LOd8ROtgpjL0XP0K1C0jjoNxSlH7eCJO92XP5fD3bwMqQh9PJMicf2YXbTMeJ3vh+Zz3V38OCkY8jhlvYAnMmdmqz5Q3ZlodaBSqetNJAxFJqepxf9eBDRv+ab7NOKr6P0Tkz6hE2L4F+HsR+RrwOVV9arrtazWs/SVwOfA9Vb1IRF5EdbbfIPYDKyc9XgEcaODxpqTo+jz6zDj/es8+XD/k5Rcs5cqNA6RikcdqRGMwHAecqYujxB0bY916Nr3rDzn6w9uwYnEWveJVmOnpc+SIYWBn52vBXwnro4RN1wERuQL4HJAGVonIhcDv9K9f3xSbvaqqiBwCDlEpc9cD/LuI3K6q7z3TtrWqlqeqR0XEEBFDVX8oIn89x3GfiXuAjSKyFngGeB3HVsBngR+ETJRc0JC4IySc2vLNjBQ8PvTNxwhCJVDlY7c/QW/a4eylaTKxDisO3mBGxwqMjhUwTYPubJJMuoHFrjsYJ52m9/Lnkb1gC4hgptItF3inQFAfIayrDsySjwMvBW4FUNWHROT5QDNs9u8B3gQMAZ8F/lhVPRExgB1AXcR+VETSVLKwfUlEjlC5qjQEVfWrixHfBUzg86q6bTb7ypd8fv7EIf7ph49R9n1+9fJ1XH/xCrLJ6UPNH9w3UhH6UMn7LqEqP37yCEVrgnP7BuhPRD7wtXBkcJwf/+wJVGFktEA6FeMV111Itqtzk7w1EjFNrExrp9CoRymHeurAHMfx9CkX1Ipdd/5n9v3Aa1R17+ROVQ1F5BXTbVyr2L8SKAF/CPwmkKUBWdkmo6rfAr411/0Mjhf54Fd/gR9U8s589NYHWDWQ4NINim2e+fZ3TX8KBcqhT1j99q7uT7J7/CimIZHY10AQVrJFveTF5xGEimEIIyN59u8/SvbcFU0eXecS+j7+xDgjd96JmUzSfcklmOnMvN0FhHUSwnrpwBx4WkSeC6iIOMB7qKxXXjDfA1HVD57hue3TbV+T2E9emABuOu0LW5C7dgweF/pj/PSxI2xZ24M9jTPCqt4Uv3ThMm6+fy8CXLlhgMvW9fDRhx5moKFua52DhkoiYTOaKx7v6+5O1hT5Oedja+VCI50TzVkz3vBRHnvfe/HGxzEQ4qtWcfb/+gusTONz4LdYLNFceQfwCSqeQfuB24B3ovx6u53kdEFVp8v8JlTWClr7XhI4e3m2ErqtJ/KRnruyC6uG3386bvGGy1fzK5csJ+95jHtFto0c5E1nbWVlJopurQXDNPDCk2vKFkouSxY37v3TMEQDj3BiiFBDJN2Hbxgk7NmvE/h+QG60gB2ziCcdzBYusBJ4Hoe+/W3Gh46iqpiGQbhvLxOPbaPnOZfPyxg6pSJntQD4b57aP7C+w4KqVLVVS+HUzLrFGd70gvP46h3b8cOAl2xZzZXnLMMyavvhp2IWKgG7cqP0OinKYwm+fO8hehJH+c3L17CiO4ldy5VjAWMaBjHHwvMDBMFxzMYWKAl9vINPUvCKKIqMHya+7GzG3QJdzszvyHJjRR57cC/nbF5CeHgfnl+AZeswkl2I1UKlAKsUPY+y6+JXTWhBtYB7OF+F3BVm4P7dkojIe1X1I6cLrupfN/9BVXOl430IMwmHN7/oLN7wgg2EGmIYIZm4g2HU/iNNOzHO6evmrp0TfPrHOwkJCDXg0YOjfPoNzyEbif1pESDmVN5r2zYRhJhjNTQeJSyMUvbLaPU3qmFAmBuhlOxiNreiT27bz+o1PQRP/BzvyH5KQHL/duJbXozRt7zlvGGGSy4911zLgdu+S1CsmM+sRYvInLd5Xo5fR2+cZvInwEeAp4ApI1brKfbV1BD3As+o6rSLrbOh48UeIJ2Y2m95JoyVivx4x0H8aiUpgJFCiacGJ7h4Ve+c99+piAiOXcno6QU+lmFiWWZjC6qLyamTMTUM3HDmDmS+53Pw6WE2ruuidGT/pP6AYP/jGNl+mIN5qBEYItw7nufSj/4doz/+EZJI0P/8F+A6Meb+S6iNei3QNpHDIrKaSvDSi6Z8RX3P8fepLPw2zDS+IMS+HjimwdLsya6ChgiLMq31Q29F/KBivhGEoudh+AGpmN0wU46RzOLEUvilicpjOw7JLIlZmBZMy6RvUeZZGRpNQ6h30el60ZtKUsbgxu27OGvzRZhicI7C2Yn5+662uxkH+BTwHWAdJ6dOEOos8yKyAng58CHgv9Vz35OJxL5G0k6CV1+0gof2j/DU4AS2afLai9fSXYe7hk4nCJWxQpFSteyiiKCaIB2PYRj1F0sFYkvOwgl9yl4JrBj5wKMvPvMlKBHhrPNWUHLLOEvX4h7cjW2bmLaNteo8sFvv8086NteevYGLVi5jtFhkcSZNX2p+vcfaXepV9ZPAJ0XkU6r6u6c+P7Bug87AjHOm/D9QCdx6L5Vswg0jEvsaSVopjGSBv37NJYzkXdIxh6RjRakTpiEIQ0TkuNBDZdZXdF3SifpHIYe+izd2kMDNYzsZYtnFhGLQX2PU9FT0LuoiN14k1nclyTWboDSBObAScRItZ68/RirmkIo5rKZ73o+t2lHeOM8S+llw2vw/1WCoI6p6n4i8sA7HOi2RUs2AuJUkbkH3PN4OdwJeEGAZxnHvEADLNCvTvzpqpR+U8UYP4JVGK4+LZZzQw+mde1bU9LFo30SUYbUW6pQbp7Wpzyk+D/hlEXkZlYzCXSLyL6r6hrrsfRKRG0lEQzENAzcISCdiVTs32KZJJhGrpVDRjDAw8MonV27yyuM1VUSKqB/HvHFqaW1NHSpVqer7VXWFqq6hkvvnB40Qeohm9hHzQHcizkSpTH+mko/IMg1MkbqbQEINMe04gVs43mdY8WokbesGQXUinWLGOSNtdo6R2Ec0HNMw6E42PulZiGJnl6EjTxP6ZcR0cHpWUHFhjphPOt2MI9Q/66Wq/gj4UX33eoJI7CM6BseKU0aJ9a8DKjlxFGlstG7E1HS21rclkdhHdBQxK0qbfCb8fK5yEQxDjFgMw65/ugelI4KqzozOfz77uRKJfUREi+KWPTzPww+FRMLGsef2c/XGx9n/hc8ycs/d2L29rPytt5DZfD5mvM7eZVop9tPxtNkpRve3EREtiJ/PkX/wZ+R/+k3Cwf2MDI6Sy5dmvb/QdTn8zf/k0I9/Qn48z+je/ez4xMcISrPf55k45ms/XWtr6uCNM59EM/uIiBYjLBU4+h//xNiuHZWO+35G36vezISxhnRqdrPwoFRi/PHtuF41P5AqxdEJigcP4nR312fgVSpmnBZSuUbRZqcYzewjIlqMIDdOYd+uEx2qFB/6OVbgnn6jaQgti8T6jSf1GYkE9sCiWe/zjMcLa2sR80dTxF5E/kZEHheRh0XkFhHprvavEZGiiDxYbf+vGeOLiGgmYpqY5sk/TbFsnNjsF1NNJ87A9a9g0QtegOE4JJYuY917/pBSWH8JUCqul7W0tqW6QFtLaxWaZca5HXh/taDwXwPvp5I/GuApVd3SpHFFRDQdI5Eife4Wgkfuxw9CxHbIXn41Zmr2NY8ty8C14sR/6XVsecvbKRXLHBkrs6oR8Q9KtEDbgjRF7FX1tkkPfwH8ajPGEdE4grAyb1OtpII2G5DdslMx4km6r/kVkpsvwx89SmztJoglseNzy7DZ25PCNA1yhTJGLMG6dT0kGpC1VVH8ts+FMD2tNGuvhVZYoH0r8K+THq8VkQeAceB/qOpPp9pIRG4AbgBYtSpKTtVKBGHISN7l6HilSlI66bAkm8CMgptqxogniK85q+77zXYlyHbNQzRzmwnhrGizc2yY2IvI94AlUzz1AVX9RvU1HwB84EvV5w4Cq1T1qIhcAvyHiJynquOn7qSaD/pGgK1bt7bZ297ZeF7IoSNjhNVfvFtyidsmPakYhkh11q9YRpTGoBPRhWLGaTMaJvaqes2ZnheRNwGvAK7WalkbVS0D5erf94nIU8BZnFwpJqLFKRbLx4UeQEOlVPLQpMMzxRESpoOJiaEGCcduSFK0diYIQ8JyEfE9DMsG08Rw2iutdjSzbz2aYsYRkeuoLMi+QFULk/oHgGFVDURkHbAR2HWa3cyY0dECh46MMZErsWggw5JFWWJz8HCImJqYY2GInORrnbBNxr0i5cBHfIvxUsXEk7JiLO5K4VjRLB9gpFQg5rmMf/9reIf2EnMSZLa+iNimSzBi7SH4Cidd7DuVdpueNMuI+vdUSnDdfoqL5fOBh0XkIeDfgXeo6nA9DpjLl8iXPQYLLhOhYtgWTz9Tl11HnIJXLLOsJ0k66ZCIWSzuTpKKWTw5foA+J814qXz8tW7oM1F0CSKna9zA5+D4CPmHf0bhwC68MKBQLpC75/uoW2z28GbEQnC9jCJoa0BVN5ym/2bg5kYcs+yH/ME/3M69TxwE4JzV/dz4Ry9jIlckk46SZ9WTZCbB8MFREqpYjoV4AYZhYBvmFJGVUi1O3W7zpPozVi4T+h7+6JHjfaGGhGFAMD6Cmelp4uhqRxW8Dp/ZNyLFcaNZMO4Rdz9+4LjQA2zfO8R37tkVeYg0ANuxWbx6gP5lPXT3Zehf1oPtWKxM9RMQkqgm9BIExzBJxmwWmmfmVHcylmFw0HWxV26a1CsYThyrd/H8DW6OKEoY1tbamjrN7EUkLiJ3i8hDIrJNRP5XI4bbCq6X88LweAnTME76kY3kSjhRwfCGYTsnr4cMxLsYcXMsyqQplj1QIRmzcUxzwSzQlstFREOCMKAUhIgdozteubPsiSdYksrgZXvouvRaik89TCKdpevSa8Guvz98I1kIuXHqOLMvAy9W1ZyI2MAdIvJtVf1F3Y7AAhL7qy9ZwydvuYeRiRIaKsmEw6uv2lQpfB3RMDTMV8sCCiI2PU6lNGHctFAqAVcLhcAt4e99lNyT96KBT2zJWszNzyfnGqSdGABn9y1iqFigfNZFDJx7KaZhYMVTbXUxPFaDtuOp0zlWvRFz1Yd2tdX9HVwwYt+TjvMvH3glX7ztEVw/4A3XnM+y/kyzh9XRhEGOQIXxoo8idCVsDA2qwiULznSjpTy57Xce/xWXD+1CehZjrDrv+Gtsw2Rpqs2/l7owZvYzoF9EJruP31iNEzqOVGpn3gdsAP5BVe+q9yAWjNjHHIsNy3t53+ufS6hKugFh4guVol/EwKLk+yRsGwEswyTE4umRoxXTmTiMlYTl3SmGckVCYHEmvqBMOGFuCu+v8aOYsxRGPywQqospcUyjddwyFfDb3R4/HZVsb7UypKpbz7g71QDYUk0KeYuIbFbVR+c0xlNYMGJ/jGQ88quvJyW/SMkL2Tk0iBeGGGKwsa+HngTky+VJayQx/BDGSx6OJQwWJgjHPZZnM9jmwvgaWj1LcewYZe+E66m1eA3WDO3xoVcmLOUp7tuG2d1L2NOD7zjErP56D3nWLISZfSOmKKo6KiI/Aq4DIrGPaB0Ug6dHJ/Cqoh5qyJ6RcbKJAQyZ9PUSA9UQESHQymsLXhkljaouiNm92DEyW6/HfvJufK9EfMU5FbG3ZjYB8QafYeiWv0erF43k+c8lecWLKYRFRCwSDagrOxMim/3MqAaTelWhTwDXAH9dn72fIBL7iDlhYFD2/ZP63DBEVUnG4jhWAtcvgYbYhkFXzGbf2AQAScch54+TtbOY0vlfRbFszP4VpLoHKh2mjczwriYs5snf84PjQo8YuNseoevil5Dffg/mig2UMlniydmnQ54zqminm3GoqzfOUuCmqt3eAL6mqv9Vt71X6fxf2BwYzZXJFT0MQ0gnbLqSkZ3/VJSQ3mSSw7nc8b6eeILKAqywoqeXguuhQNK2GSkWCVTJxGL0JG1GvCP0xvqaNv75RgwDnLkF8anrHf/bIkFYKhDkJ5i4/w7yt/0Hy3/jd9FVG5p2t7RgZvZ1QlUfBi5q9HEisT8NIxMlPnDTXfzgwf2YJvzyFet4/2sviQT/FOJWjJXdJrZpMF7ySDsOy7qS2NVKS6YImXjFrdANyjiOx3Inwbg/yr7iEGtTG8+0+4hTkFiC5PmXUz64A1UfFJxla0CE4uAgGgZM3P0j4gPLsJo4u18INvtWSoVQC1H46BR4fsDt9z/N9x94GlXF95Vb7niKh3YNNXtoLYljWizvSnPWQDcru9OnTWrmmDHSVgaVgLSdZlNmM0krPc+jbW/EMIitPYe+V72D9HlXkb7yZWRf+noGf/CfBGEAVGIXmrkGogpeENbU5oKI/Fo14jQUka2nPPd+EdkpIk+IyEsn9V8iIo9Un/u/Mts3SkHC2lqrEM3sp6BU9tl16OQU+qrw1MExrtq8rEmjam1Mw6CW8DTTMOkyuhs9nI7GiMWJrToLZ+kawiBg6Dv/Rn7fThCwLYfspS/ATCSbNr6KGWdepr2PAq8BPj25U0TOBV4HnAcsA74nImdV3Rs/RaXo0S+Ab1Hxevn2fAy22URiPwWWIVy9ZTk33f748dtR2zR4wfnLG3rcsudS9pWxkkd3wiFuW8fNIRERpyK2g2lD/0t/hfSGc/BHhkiccxFWV7MTpinBPCzQqup2YKq7mFcCX63Wx9gtIjuBy0RkD9ClqndWt/si8CpmIfbtmAgtEvspSCQc1izK8De//Vy+/KMnsU2Dt770XPozsYYd0/VKPHpgjE/+cBu5skt3Mskfv+RCNg50YSy0UNMm4fo+oor6LhL6mPEUIkZlUbWFMRMpUpsvbfYwjqM6o3z200aXzoLlVGbux9hf7fOqf5/aPzsise8M+rJJnn/eErZuHMAwhGzSaWihEzeAz9yxnVzZBWC0UOALP9/Bn15/AV1zLDQdMT1hqHiuRzgxiJsbQ4CY45BYshqZh+hU9T3QECynI2IOZjCzP2N0aS3lTafabIq+0+XRbjPJnj2R2J8GwxC6Mgm65uuAAkdzpZO6Do3lO2oF3Q89glCrNl0PUwziVmvUEhgvFLEIcXNjQEUByq5LbGIY6V6MNKheroYhWsrj734IdfNYS9di9PaCKYjZni6plQXa+mjodOVNT8N+YOWkxyuAA9X+FVP0z2JgIG0WS9BJWtLWqAZsXn7yj/uiVQMdMcuDSv72oheya2ScJwaHGcy5KIIbuM0eGlCpmRoEJweHKRUxbiiBi3vftwmeeYzwyHbch75FOHwIde9Dg5HGHrtBKJWZfS2tQdwKvE5EYiKylkp507tV9SAwISKXV71wfgs43d1BxxHN7FuEtGPy7hedw5fv2cPuwXHOXdrLr16yllSH1MgNVdk9MoZXFdSjhSKGCAPpxq2DzISYbRIGDqYTJ3Ard1imIRipLEjj5kThxChaylFJaV7BP/AUdmYxhEfAbPZi68xRnZ8FWhF5NfBJYAD4pog8qKovVdVtIvI14DHAB95Z9cQB+F3gn4AElYXZWXvitJJbZS1EYt8iiBGnJ+nzlivW4oeCYxoknNYQwnrghyF+GJzUl3c9FlEx43hBUH1NSKDh8YIe80UqHmNkokC8fyVhYRRDA+xUFiOWaOjdlUxRRFxiSaAA2sSUB3NknrxxbgFuOc1zHwI+NEX/vcDm+gygLnuZN5pixhGR/ykiz1SLjT8oIi+b9NyUwRALARGLVCxFNpEk4bROytp6YBsmtnHy3CJpVxYj/TBkuFBm70iOoUIZEeFgbvw0e2ocPZkkhmkiqR6s7CKsRLphtvpjiJPAXH42UFmEl3gX1przQJ8Cc2lDj90oWsCMMy9Ija1VaObM/u9U9W8nd0wTDLEgccsehVyZRCpGrI3TMxuGsLo7y/6xHOXAoysWY2kmhSkmQ/kS+8fzAOQ9H9cPySZNcp5Lep7L8cWc+X2PxY5hb7wEa/W5aHkCI5VG9UnEfjlitk7K4hmh8xZU1Vza7BxbzYwzZTAEcGdzh9UcBg+O8tPvPsqh/cMsWd7DlS89n0XLups9rFlhiJB2bM4a6K50qGKZJkGoTJRPXqTNuR7LsgnKvjfvYt8MxI4hdgxS3WhYRmQrlQSI7Ymic06F0PLUWEy8lWim2L9LRH4LuBf4I1Ud4fTBEM9CRG6gEvbMqlWrGjzU+Sc3XuRbX7ubctHjeS+/kCUreimWPfL5EqlUe5p4DEMwnnVjqzimheAe/+1YhiAIyRkKvRd4hATk/AlMMclYWcwGm2HqjRgzX6dRz0UnrYcYsea6s6rOj82+2UQRtFXOFAxBJT/FX1C5Nv4F8FHgrcwg6KEaZXcjwNatW9vsbZ+eUtHFLbm87M1X8YW793H/T/ewsjfFu6/dxNq4g9UhaRSOeeSUfJ+85yEoyzLJSulIqzaxV1VCVYquixt42FYMFZ+DxadZkVrT2BNoMqFbwn3qUfL3/wDcMs6qs0k97+UY8SbnxlkAYh/N7KvUGgwhIp8BjiXqP10wxILDMAy2Xn0eX7n/Ge54chCAJw9P8Lffe5K//ZULyHSI2IsIMdNkfV+GchBiGYKqErNq/2oqkCsWKfgFALwAkrE4OX+cgp/r6MyaQanA0Z/dShD6CEK4ZxtW/xLim6+YcWGUetLpWl/PxVcRWQl8kcrkOKSSMuITddr9cZrljTPZzeDVnKi1OGUwxHyPrxXo6kmybO0iHjtwzCtFcGI2Q3mXotdZ69WGIZiGQdK2cExzRkIPVb9uPdlG7Pk+fbFF+OqfZqv2pxS45I8+UzFfVd+DYuDiDT5TSb/QJCpmnLCm1tZojW16fCqm7HOAy4F3Vp1V6kqzLv0fEZEtVN6KPcDvAEwTDLGgsCyTRNLm3NU9DLtBxd5tCL0ph4TdXnboRiNUyiOKVO4KAEwxCBCSZvv6qk9H3ivj9C1FTAutBquFqljL11cWfJtG+7tV1kSdTrEa2Xuw+veEiGynslb5WH2OUKEpYq+qbzzDc1MGQyxEuhIOb7liDblywGMHx1ncFecdV60lfpriIAsWERKOg7qKF7qICIlYjBATy2hfd9XpsAyTPeUxVlz9Wgr3/oCwlCe5cQv22nObmqlTFVy/zWft06E6E1tVzZk9RWQNlRKFd81tgM+m1VwvI06hL+Xw36/ZiCFCqErCNjtmcbZeGCLYtkXWMgn1WMSrYs5DtspmknWS2E6cu51xLrj2tSRNm6IGhE12V53H4iVNZQbeOGfM7Hl8fyJp4GbgD1S17lGFkdi3OCJCpo2DqeYLQwREFlxmv42ZJfTFMkx4RcSy6HWyxMzmf19mkM++janfOYqITUXov6SqX6/bjicRiX1ERBtjGiYD8S4G4vOWjHta5isRWrOpl599NQPn54Dtqvqx+uz12URiHxERUVciP/sZ8zzgjcAjIvJgte9PVfVbdTsCkdhHRETUG4WgTsVLWhalbsEEqnoH85AzLRL7iIiIuqIsDJt9lC4hIiJiYaMcj3eIaB0isV8AuH5IwQt44sgE/akYS7tiJJ3oo49oFIougJl9lOI4ouU4MF7kH+/YQ9mvBCNftqqHV1+wNBL8iIYQmXFak4XmlrzgKLg+tz1+hILrUfA8Cq7HnXuOUuiw/DoRLYRWJr21tPamfslx5oNoatfhhApFPyDvemj1i+e7IaVI7CMaSNjhxUtE26/geDSz73CStsnlq3uYPMNY35/CMjRaRJtnVENU20whZoOChlpTa2vaa2Ifzew7HcMQlmcd3n3Veh48MM6STIxLV3YzlM/Rk7KwDRtDomt+I1ENUd9DJ4ZAFaNrAEy7qcnKGokuiAXa9rNDRWK/ADAFHj/8NJetXMrSbJYjuTI9qRSHC0Mc8fawIX0O3U5vs4fZsYS+T2nPo/heuVKsZewIzsrzYBYlCNuFjtd6iMQ+ovVYms1w/bmbGC8FfOT7TzJaKhHi8fJzl3P2in5+OngbL178clJWptlD7Tg838cfH8Itl6o9SrFYwpoYwuhZVs3Q2Xl0vImwxUw0tdCZ95ERz8KxbL77+CA5NyDAJ1T41vYD9NhLGPWGKQaFZg+xIyl5AaeuVaoq2vjo+OaxQGz2orW1ViGa2S8QwlAZL/lViankew9CJe/6JMwUhkQFURqBAUiqG8NJELrFSp/jIJm+zp3V0/neOAC02WJ7JPYLhIRtcsnKbvaO5LHFxtUyK7qTJB1YEl/W0eX76kUQhuTLLinHQlQrC6wiyBkWuFOJGIeOjpNZfjZBfhQBYl29BGrSsZdXbTsdnB0tNGuvhaaIvYj8K7Cp+rAbGFXVLdWSXNuBJ6rP/UJV3zH/I+w8LNPgOat7iFkGDzwzSl/S5upNfYTGOBf3XEHcTDR7iC1N0fV48sgQ6/uyjB09jEFI3LKxU10YTvyMgr+4N8NIrkAQ68I2TULDJN7R0cu6QCJo2+scm1WD9teP/S0iHwXGJj39lKpumfdBLQCSjsVzVvdwwbIstiHEbBOIZvS1cGhiAlMgKE4Q+B4BlepYRmECw46dMUGtiNCbWTjvsxIt0LYiTV2grVZoeS3wlWaOYyFhGgbpmFUV+oha2XN0lK54jMD3j/cFYYgStp0LXsNZEAu09csJISKfF5EjIvJoI0fcbG+cq4DDqrpjUt9aEXlARH4sIlc1a2AREZNZls2wY/AoVuyEucsQQQwLoqC0ZxEEYU2trQm1tjY9/wRc19jBNtCMIyLfA5ZM8dQHVPUb1b9fz8mz+oPAKlU9KiKXAP8hIudNVWldRG4AbgBYtWpVfQffYLwgpFDyiTtmNMNuE5Zlu9g5OIxrOMRSXeCXicfi2IkMdKhXzWzRhbJAWydU9SfV9cqG0jCxV9VrzvS8iFjAa4BLJm1TBsrVv+8TkaeAs4B7p9j/jcCNAFu3bm2L+0E39MkXff7z7n3c8+QRNq/q5bVXracn3bmRlJ1CJh7jxWetY7xcBidJItWFZRgdm/JgbiyQvEttdo7NdAm4BnhcVfcf6xCRAWBYVQMRWQdsBHY1a4D1xA19Dk6M8c+37+aff/QEhgh3bD/EjgPjfPB1F9OVcpo9xIhpSDg2Ccdu9jDago4X+8oqdK2v7heRyRPWG6uT1XmlmWL/Op69MPt84H+LiA8EwDtUdXjeR9YARt08hpp8+769AISqlAOPnz52AOTiJo8uIqKOVBdoO57aT3FIVbc2cCQ10TSxV9U3T9F3M3Dz/I+m8ZQDj4If0tcVZ2iikiclRMkkoxl9RGexMCpVtV/Wy8jgOE8krBi7Sgd41ys2Hy8HGDNN3v3yzdhW9DFEdBCqhL5XU2trwrC2Ng0i8hXgTmCTiOwXkbc1YridHMbXUvQ6KdZ3LWbYLvCN/3EdOw+OcfbSPrqSDomOjqaMWHgoqh1eCW1mNvsz70r19XXZ0TREKlMDqkqx7OJYYFmz85wxxOCsrqWMxQv4YcAVm5aQsM5swim6PiUvxDQEMX3yQRlLDLqdFJYRuWxGtCgKGjRe7EXkb4BfAlzgKeAtqjpafe79wNuorP29R1W/W+2/hIpfewL4FvD7OuvV5PYy40RiPw1juQI/uXc7P7xnG1vPW8V1V15AOmXjmLPL/Z51krUdt+DyLz/ZxZ1PDrKsN8Fbr1nHfYXH2Jk7yHXLtnBR71ocI/r4IloRrcl8UQduB96vqr6I/DXwfuBPRORcKg4g5wHLgO+JyFlaud34FJX4nF9QEfvrgG/P6uhtFksQqcUZyBfLfObff8Tf/ct/oOrz1e/A/duv4r1veyH9XXEMozFueCXX58s/3c1Xf7YHBbYfHObxQyN88rcv5cHRPfzX/vtZk17E4ni2IcePiJgLCmjY+Jm9qt426eEvgF+t/v1K4KvVuJ3dIrITuExE9gBdqnongIh8EXgVsxL79kuOE60MngHPD/jqt3+O6ol8KDfffie2mcUNx86w5dwoeQG/2DEIVExIIcreoRy5ok+3nSQflMl7pWn2EhHRJLRis6+lUfVBn9RumOVR38oJ0V4OPD3puf3VvuXVv0/tnx11yo0zX0Qz+zOgqiTiMY4V+wCIxxwqdYYad500DYOVfSl2Hc4hIhgYJGJCNumQ98tkrAQZO0pJ3GkErkvolgmDEDMeB8PEstvxJ6oz8bQ5ow96LWlXROQDgA986dhmUw7q9P0zp+JfOqtNm0U7fpPmjWTc4fff8FL+8G++QBCUAXj3b1yHkscxVzbsuJmEzTtesol9Q3l2H8nRn0rwe9dvZGf+AFknyXXLttDjNDZlbsH1KfshfqgkbZNUzKKULzExnCfwAxKZBNn+Z69bBGFIqGCbBhqWQSwkqoI1LaHn4R7YSXlkBFCsdBfJNWcTlEuIBogdA8Nsi+pWqkpYJzNODWlX3gS8Arh60kLrfmDyD3QFcKDav2KK/tkObtabNoNI7M9AzLF56XMv4Eef+5/c9cgOLty0nGWLsmRSiYb/6Jb2xPn7t11GruyTdCzECClLmQv6VtDjpBp6/ILr8/WHDnDnnmEU2LQozVsvW81937yPn958F17JZeXZK3jlO69j0ap+oPIDL3gBjxwcZ7zkcd7iBP3xo9jhQxC7HDGnmpxFHCMsjFeFvpJEzMaD/Ah+fgJ/5BCxTAZz+TkwjQdXqzAfNnsRuQ74E+AFqjq5iPKtwJdF5GNUFmg3AndX07BMiMjlwF3AbwGfnPUA2kvrI7GfjkwqTiYVZ/3KxfN6XNMw6Eo6dJ0UYRufl2MfGCvx8z0nslQ8cSTHXXuGKY0VcYsuAPu27+fn/3kPL3/7NdiOTdkPufnhAwzlXaDMtoMlXrV5Bavjw2j+Jki9HTF752X8bckxs4cqYgimaaBeGTFNPC/AKRXQkQPQtxJpebdbRecn7eXfAzHg9urk5xeq+g5V3SYiXwMeo2LeeaeecPz/XU64Xn6b2XritGEEbST2Ec/i0ET5WX0Hxoqs7D7ZdDS47yjlQhnbsRktelWhD0FdFHjo4ARL1m8h5n0awmGIxP60GOluDNskdH3EqJaFT3XjHdyFohXfD6/cHqXw5snPXlU3nOG5DwEfmqL/XmBznQZQl93MF5HYRzyLcxdnsA3Bm5Tf5NLVPTz0owdOet3qc1eQSFcWim3zmFlJOLYOZpsGggsYINFX7UyIbZNatxlv+CAahBiLlhGUCrjjo5imgYgg6T5o+Vk9FW+cdk+FMB3RAm1EJ5CKmbzzqvV8Z/sh3EB5/vp+VvUkyW1Zw+G9gxTzJTZetJZLr9uCaVXEJ+VYnDWQ5snBHIhDzHS5bGUXTnAr2BeC0dPks2ptxDAxk2mw1+B7AcMFl6QVJ9G3GIsAIzuApLrbY4F2IaRLgGhmH9H+xCyTjQMplmXXAJCwDUzD4OJrL2DDRWsJg5BUNkk8dWINIW6bXLNxgC3LskyUPVZ1x3B0F2I8F8yNiBEFgNWCaduYto0Vc3A9H2fZhsp9kmkibVT+UNts1jtzIpt9RIcgIqRjJ389TNOkZ3H3abeJ2yYruhNU1r4AzmnY+Dod0zRImO3hefMsVOfFG6fZtFuBlkjsIyIi6k7Hz+zrmPVyvojEPiIioq7UM6iqpYnEPiIiYkGjM0qX0J6ookF73b20z4pPREREm6AQBrW1tkZrbGdGRK4TkSdEZKeIvK9Ro41m9hEREXVnISzQUoc6u1JJHPUPwLVUcvfcIyK3qupjc975KURiHxERUV8WiDdOnZLjXAbsVNVdACLyVSr5+Osu9tJu7kNTISKDwN5mj6PO9ANDzR5Eg+jUc+vE81qtqgMz2UBEvkPlvaiFIVW9bubDai4zPMc4MLkAxY2qemN1P78KXKeqv119/EbgOar6rnqOFzpkZj/TL2M7ICL3ninPdzvTqefWqec1U9pRvGdKHc+xfjn2pyFaoI2IiIhoHqfLvV93IrGPiIiIaB73ABtFZK2IOFQKpd/aiAN1hBmnQ7mx2QNoIJ16bp16XhENQlV9EXkX8F3ABD6vqtsacayOWKCNiIiIiDgzkRknIiIiYgEQiX1ERETEAiAS+yYjIr8mIttEJBSRrac89/5qCPUTIvLSSf2XiMgj1ef+r7RDRQvmLyy8UYjI50XkiIg8OqmvV0RuF5Ed1f97Jj035ecXEdEMIrFvPo8CrwF+MrlTRM6lsjJ/HnAd8I/V0GqATwE3ABurreX9mieFhV8PnAu8vnqO7cQ/8ez3+n3A91V1I/D96uPpPr+IiHknEvsmo6rbVfWJKZ56JfBVVS2r6m5gJ3CZiCwFulT1Tq2srn8ReNX8jXjWHA8LV1UXOBYW3jao6k+A4VO6XwncVP37Jk58FlN+fvMxzoiIqYjEvnVZDjw96fH+at/y6t+n9rc6pzufdmexqh4EqP6/qNrfqecb0aZEfvbzgIh8D1gyxVMfUNVvnG6zKfr0DP2tTruOe7YstPONaHEisZ8HVPWaWWx2ujDq/dW/T+1vdeYtLHyeOSwiS1X1YNXEdqTa36nnG9GmRGac1uVW4HUiEhORtVQWYu+umgomROTyqhfObwGnuztoJeYtLHyeuRV4U/XvN3His5jy82vC+CIigGhm33RE5NXAJ4EB4Jsi8qCqvlRVt4nI16jktfaBd6rqsSThv0vFMyQBfLvaWpr5DAtvFCLyFeCFQL+I7Af+HPgw8DUReRuwD/g1gGk+v4iIeSdKlxARERGxAIjMOBERERELgEjsIyIiIhYAkdhHRERELAAisY+IiIhYAERiHxEREbEAiMQ+oiGISK7ZY4iIiDhBJPYRERERC4BI7CMailT4GxF5tJqD/9er/S8UkR+JyL+LyOMi8qV2ycsfEdGORBG0EY3mNcAW4EKgH7hHRI7l7r+ISr73A8DPgOcBdzRhjBERHU80s49oNFcCX1HVQFUPAz8GLq0+d7eq7lfVEHgQWNOcIUZEdD6R2Ec0mjOZZsqT/g6I7jQjIhpGJPYRjeYnwK+LiCkiA8DzibI/RkTMO9FMKqLR3AJcATxEpXjHe1X1kIic3dxhRUQsLKKslxERERELgMiMExEREbEAiMQ+IiIiYgEQiX1ERETEAiAS+4iIiIgFQCT2EREREQuASOwjIiIiFgCR2EdEREQsAP4/qaBOkhii+sYAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], "source": [ "ds_mesh2.plot.scatter(x='longitude', y='latitude', c='k', alpha=0.7);\n", "ds_selection.plot.scatter(x='lon', y='lat', hue='field', alpha=0.9);" ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { "kernelspec": { - "display_name": "Python [conda env:xoak_dev]", + "display_name": "Python 3.9.13 ('extract_model')", "language": "python", - "name": "conda-env-xoak_dev-py" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -195,7 +374,12 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.6" + "version": "3.9.13" + }, + "vscode": { + "interpreter": { + "hash": "8838a114cb9af02008a6426a5905618667030dfe4375f4e8b79cae9729b737db" + } } }, "nbformat": 4, diff --git a/src/xoak/accessor.py b/src/xoak/accessor.py index 12e52f4..95639d8 100644 --- a/src/xoak/accessor.py +++ b/src/xoak/accessor.py @@ -1,4 +1,4 @@ -from typing import Any, Hashable, Iterable, List, Mapping, Tuple, Type, Union +from typing import Any, Hashable, Iterable, List, Mapping, Optional, Tuple, Type, Union import numpy as np import xarray as xr @@ -134,12 +134,18 @@ def index(self) -> Union[None, Index, Iterable[Index]]: return [wrp.index for wrp in index_wrappers] def _query(self, indexers): + """Find the distance(s) and indices of nearest point(s). + + Note that the distance is converted in function from radians to kilometers + by multiplying by the radius of the earth and dividing the resulting meters by 1000. + """ X = coords_to_point_array([indexers[c] for c in self._index_coords]) if isinstance(X, np.ndarray) and isinstance(self._index, XoakIndexWrapper): # directly call index wrapper's query method res = self._index.query(X) results = res['indices'][:, 0] + distances = res['distances'][:, 0] else: # Two-stage lazy query with dask @@ -195,7 +201,7 @@ def _query(self, indexers): concatenate=True, ) - return results + return results, distances*6371000/1000 def _get_pos_indexers(self, indices, indexers): """Returns positional indexers based on the query results and the @@ -225,7 +231,7 @@ def _get_pos_indexers(self, indices, indexers): return pos_indexers def sel( - self, indexers: Mapping[Hashable, Any] = None, **indexers_kwargs: Any + self, indexers: Mapping[Hashable, Any] = None, distances_name: Optional[str] = None, **indexers_kwargs: Any, ) -> Union[xr.Dataset, xr.DataArray]: """Selection based on a ball tree index. @@ -242,7 +248,19 @@ def sel( This triggers :func:`dask.compute` if the given indexers and/or the index coordinates are chunked. - + + Parameters + ---------- + distances_name: str, optional + If a string is input, it is used to save the distances into the xarray + object. Distances are in km. + + Returns + ------- + xr.Dataset, xr.DataArray + Normally, the type input is the type output. However, if you input a + str for `distances_name`, the return type will be Dataset to + accommodate the additional variable. """ if not getattr(self, '_index', False): raise ValueError( @@ -250,7 +268,7 @@ def sel( ) indexers = either_dict_or_kwargs(indexers, indexers_kwargs, 'xoak.sel') - indices = self._query(indexers) + indices, distances = self._query(indexers) if not isinstance(indices, np.ndarray): # TODO: remove (see todo below) @@ -262,5 +280,14 @@ def sel( # as OuterIndexer, while we want here VectorizedIndexer # This would also allow lazy selection result = self._xarray_obj.isel(indexers=pos_indexers) + + # save distances as a new variable in xarray object if name is input + if distances_name is not None: + # need to have a Dataset instead of DataArray to add a new variable + # otherwise goes in as a coordinate + if not isinstance(result, xr.Dataset): + result = result.to_dataset() + # use same dimensions as indexers + result[distances_name] = (indexers[list(indexers.keys())[0]].dims, distances) return result diff --git a/src/xoak/tests/conftest.py b/src/xoak/tests/conftest.py index e8b4314..3774cba 100644 --- a/src/xoak/tests/conftest.py +++ b/src/xoak/tests/conftest.py @@ -1,4 +1,4 @@ -import dask +import dask.array import numpy as np import pytest import xarray as xr diff --git a/src/xoak/tests/test_accessor.py b/src/xoak/tests/test_accessor.py index 32877ac..651666e 100644 --- a/src/xoak/tests/test_accessor.py +++ b/src/xoak/tests/test_accessor.py @@ -1,4 +1,5 @@ import pytest +import numpy as np import xarray as xr from scipy.spatial import cKDTree @@ -80,3 +81,23 @@ def test_index_property(): ds_chunk = ds.chunk(2) ds_chunk.xoak.set_index(['x', 'y'], 'scipy_kdtree') assert isinstance(ds_chunk.xoak.index, list) + + +def test_distances(): + + ds = xr.Dataset( + coords={ + 'x': ('a', [0, 1, 2, 3]), + 'y': ('a', [0, 1, 2, 3]), + } + ) + + ds_to_find = xr.Dataset({"lat_to_find": ("a", [0, 0]), "lon_to_find": ("a", [0, 0.5])}) + ds.xoak.set_index(['y', 'x'], "sklearn_geo_balltree") + + output = ds.xoak.sel( + {'y': ds_to_find.lat_to_find, 'x': ds_to_find.lon_to_find}, distances_name="distances" + ) + + assert isinstance(output, xr.Dataset) + assert np.allclose(output["distances"], [0, 55.59746332]) From 4f3ecc64ba8a73b6280ba8ef50b33388d66be692 Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 19 Jan 2023 12:40:01 -0600 Subject: [PATCH 02/14] got precommit needed some updates in precommit too --- .pre-commit-config.yaml | 4 ++-- src/xoak/accessor.py | 17 ++++++++++------- src/xoak/tests/test_accessor.py | 14 +++++++------- 3 files changed, 19 insertions(+), 16 deletions(-) diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 284a0ae..cb23a2b 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -9,12 +9,12 @@ repos: - id: double-quote-string-fixer - repo: https://github.com/ambv/black - rev: 20.8b1 + rev: 22.3.0 # https://github.com/psf/black/issues/2964 hooks: - id: black args: ["--line-length", "100", "--skip-string-normalization"] - - repo: https://gitlab.com/PyCQA/flake8 + - repo: https://github.com/pycqa/flake8 rev: 3.8.4 hooks: - id: flake8 diff --git a/src/xoak/accessor.py b/src/xoak/accessor.py index 95639d8..5f40153 100644 --- a/src/xoak/accessor.py +++ b/src/xoak/accessor.py @@ -135,8 +135,8 @@ def index(self) -> Union[None, Index, Iterable[Index]]: def _query(self, indexers): """Find the distance(s) and indices of nearest point(s). - - Note that the distance is converted in function from radians to kilometers + + Note that the distance is converted in function from radians to kilometers by multiplying by the radius of the earth and dividing the resulting meters by 1000. """ X = coords_to_point_array([indexers[c] for c in self._index_coords]) @@ -201,7 +201,7 @@ def _query(self, indexers): concatenate=True, ) - return results, distances*6371000/1000 + return results, distances * 6371000 / 1000 def _get_pos_indexers(self, indices, indexers): """Returns positional indexers based on the query results and the @@ -231,7 +231,10 @@ def _get_pos_indexers(self, indices, indexers): return pos_indexers def sel( - self, indexers: Mapping[Hashable, Any] = None, distances_name: Optional[str] = None, **indexers_kwargs: Any, + self, + indexers: Mapping[Hashable, Any] = None, + distances_name: Optional[str] = None, + **indexers_kwargs: Any, ) -> Union[xr.Dataset, xr.DataArray]: """Selection based on a ball tree index. @@ -248,7 +251,7 @@ def sel( This triggers :func:`dask.compute` if the given indexers and/or the index coordinates are chunked. - + Parameters ---------- distances_name: str, optional @@ -260,7 +263,7 @@ def sel( xr.Dataset, xr.DataArray Normally, the type input is the type output. However, if you input a str for `distances_name`, the return type will be Dataset to - accommodate the additional variable. + accommodate the additional variable. """ if not getattr(self, '_index', False): raise ValueError( @@ -280,7 +283,7 @@ def sel( # as OuterIndexer, while we want here VectorizedIndexer # This would also allow lazy selection result = self._xarray_obj.isel(indexers=pos_indexers) - + # save distances as a new variable in xarray object if name is input if distances_name is not None: # need to have a Dataset instead of DataArray to add a new variable diff --git a/src/xoak/tests/test_accessor.py b/src/xoak/tests/test_accessor.py index 651666e..6492c7e 100644 --- a/src/xoak/tests/test_accessor.py +++ b/src/xoak/tests/test_accessor.py @@ -1,5 +1,5 @@ -import pytest import numpy as np +import pytest import xarray as xr from scipy.spatial import cKDTree @@ -92,12 +92,12 @@ def test_distances(): } ) - ds_to_find = xr.Dataset({"lat_to_find": ("a", [0, 0]), "lon_to_find": ("a", [0, 0.5])}) - ds.xoak.set_index(['y', 'x'], "sklearn_geo_balltree") - + ds_to_find = xr.Dataset({'lat_to_find': ('a', [0, 0]), 'lon_to_find': ('a', [0, 0.5])}) + ds.xoak.set_index(['y', 'x'], 'sklearn_geo_balltree') + output = ds.xoak.sel( - {'y': ds_to_find.lat_to_find, 'x': ds_to_find.lon_to_find}, distances_name="distances" + {'y': ds_to_find.lat_to_find, 'x': ds_to_find.lon_to_find}, distances_name='distances' ) - + assert isinstance(output, xr.Dataset) - assert np.allclose(output["distances"], [0, 55.59746332]) + assert np.allclose(output['distances'], [0, 55.59746332]) From 88fadb7eb5d255aa4df6e44366372fa588f62ba3 Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 19 Jan 2023 12:48:20 -0600 Subject: [PATCH 03/14] realized I should add some attributes --- src/xoak/accessor.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/xoak/accessor.py b/src/xoak/accessor.py index 5f40153..c876fcf 100644 --- a/src/xoak/accessor.py +++ b/src/xoak/accessor.py @@ -292,5 +292,8 @@ def sel( result = result.to_dataset() # use same dimensions as indexers result[distances_name] = (indexers[list(indexers.keys())[0]].dims, distances) + # add attributes + result[distances_name].attrs["units"] = "km" + result[distances_name].attrs["long_name"] = "Distance from location to nearest comparison point." return result From e7e642906f9518c873aaf4542b75b561bacf730f Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 19 Jan 2023 12:50:01 -0600 Subject: [PATCH 04/14] recently have been issues with mamba https://github.com/conda-incubator/setup-miniconda/issues/274#issuecomment-1384764160 --- .github/workflows/test.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/test.yaml b/.github/workflows/test.yaml index d0a1686..2b43c56 100644 --- a/.github/workflows/test.yaml +++ b/.github/workflows/test.yaml @@ -24,7 +24,7 @@ jobs: uses: conda-incubator/setup-miniconda@v2 with: python-version: ${{ matrix.python-version }} - mamba-version: "*" + miniforge-variant: Mambaforge channels: conda-forge,defaults channel-priority: true auto-activate-base: false From 9e693dec51f447b6f2495d65480564b3b636d90f Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 19 Jan 2023 12:52:50 -0600 Subject: [PATCH 05/14] had more lint from changes --- src/xoak/accessor.py | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/xoak/accessor.py b/src/xoak/accessor.py index c876fcf..bfb8001 100644 --- a/src/xoak/accessor.py +++ b/src/xoak/accessor.py @@ -293,7 +293,9 @@ def sel( # use same dimensions as indexers result[distances_name] = (indexers[list(indexers.keys())[0]].dims, distances) # add attributes - result[distances_name].attrs["units"] = "km" - result[distances_name].attrs["long_name"] = "Distance from location to nearest comparison point." + result[distances_name].attrs['units'] = 'km' + result[distances_name].attrs[ + 'long_name' + ] = 'Distance from location to nearest comparison point.' return result From 77acd76849bd0080faf6dbbe2fe175990d006802 Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 19 Jan 2023 13:11:03 -0600 Subject: [PATCH 06/14] trying building docs with mamba --- .readthedocs.yml | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/.readthedocs.yml b/.readthedocs.yml index 45f714a..aaa4260 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -13,3 +13,8 @@ python: install: - method: pip path: . + +build: + os: "ubuntu-20.04" + tools: + python: "mambaforge-4.10" From 524d791e79e4ed41cc97d8972b5988a349dda059 Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 19 Jan 2023 13:12:49 -0600 Subject: [PATCH 07/14] another readthedocs change --- .readthedocs.yml | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/.readthedocs.yml b/.readthedocs.yml index aaa4260..351eed3 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -1,6 +1,11 @@ # Required version: 2 +build: + os: "ubuntu-20.04" + tools: + python: "mambaforge-4.10" + # Build documentation in the doc/ directory with Sphinx sphinx: configuration: doc/conf.py @@ -9,12 +14,6 @@ conda: environment: environment_doc.yml python: - version: 3.8 install: - method: pip path: . - -build: - os: "ubuntu-20.04" - tools: - python: "mambaforge-4.10" From 7ce86c654862b02bf0a540e162fd649264084763 Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 19 Jan 2023 14:38:08 -0600 Subject: [PATCH 08/14] rearranged and matched distances shape to dims --- src/xoak/accessor.py | 9 +++------ 1 file changed, 3 insertions(+), 6 deletions(-) diff --git a/src/xoak/accessor.py b/src/xoak/accessor.py index bfb8001..6ec0b0b 100644 --- a/src/xoak/accessor.py +++ b/src/xoak/accessor.py @@ -291,11 +291,8 @@ def sel( if not isinstance(result, xr.Dataset): result = result.to_dataset() # use same dimensions as indexers - result[distances_name] = (indexers[list(indexers.keys())[0]].dims, distances) - # add attributes - result[distances_name].attrs['units'] = 'km' - result[distances_name].attrs[ - 'long_name' - ] = 'Distance from location to nearest comparison point.' + attrs = {'units': 'km', 'long_name': 'Distance from location to nearest comparison point.'} + dims = indexers[list(indexers.keys())[0]].dims + result[distances_name] = (dims, distances.reshape(len(dims)), attrs) return result From 70d00b3f0f6058c9d374dec3040ab9fdf7ed0bab Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 19 Jan 2023 14:47:45 -0600 Subject: [PATCH 09/14] forgot lint again --- src/xoak/accessor.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/xoak/accessor.py b/src/xoak/accessor.py index 6ec0b0b..9143756 100644 --- a/src/xoak/accessor.py +++ b/src/xoak/accessor.py @@ -291,7 +291,10 @@ def sel( if not isinstance(result, xr.Dataset): result = result.to_dataset() # use same dimensions as indexers - attrs = {'units': 'km', 'long_name': 'Distance from location to nearest comparison point.'} + attrs = { + 'units': 'km', + 'long_name': 'Distance from location to nearest comparison point.', + } dims = indexers[list(indexers.keys())[0]].dims result[distances_name] = (dims, distances.reshape(len(dims)), attrs) From f8bc24647d81dd8d4b3a23d6ee17b871f24c028f Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 19 Jan 2023 15:01:30 -0600 Subject: [PATCH 10/14] found a nuance while using code --- src/xoak/accessor.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/xoak/accessor.py b/src/xoak/accessor.py index 9143756..9c1e770 100644 --- a/src/xoak/accessor.py +++ b/src/xoak/accessor.py @@ -296,6 +296,9 @@ def sel( 'long_name': 'Distance from location to nearest comparison point.', } dims = indexers[list(indexers.keys())[0]].dims - result[distances_name] = (dims, distances.reshape(len(dims)), attrs) + if distances.ndim != len(dims): + result[distances_name] = (dims, distances.reshape(len(dims)), attrs) + else: + result[distances_name] = (dims, distances, attrs) return result From 611607ccb5dae82fda737f609f2b383539048a78 Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 19 Jan 2023 15:27:49 -0600 Subject: [PATCH 11/14] figured out a better way to shape distances from the code --- src/xoak/accessor.py | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/src/xoak/accessor.py b/src/xoak/accessor.py index 9c1e770..201567c 100644 --- a/src/xoak/accessor.py +++ b/src/xoak/accessor.py @@ -295,10 +295,11 @@ def sel( 'units': 'km', 'long_name': 'Distance from location to nearest comparison point.', } - dims = indexers[list(indexers.keys())[0]].dims - if distances.ndim != len(dims): - result[distances_name] = (dims, distances.reshape(len(dims)), attrs) - else: - result[distances_name] = (dims, distances, attrs) + + indexer_dim = list(indexers.values())[0].dims + indexer_shape = indexers[list(indexers.keys())[0]].shape + result[distances_name] = xr.Variable( + indexer_dim, distances.reshape(indexer_shape), attrs + ) return result From c8c34117af0ddcd0fd6fb22e46cb02a79fd69084 Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Thu, 19 Jan 2023 15:42:19 -0600 Subject: [PATCH 12/14] changed multiplication to km not meters divided by 1000 --- src/xoak/accessor.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/xoak/accessor.py b/src/xoak/accessor.py index 201567c..f7c1294 100644 --- a/src/xoak/accessor.py +++ b/src/xoak/accessor.py @@ -137,7 +137,7 @@ def _query(self, indexers): """Find the distance(s) and indices of nearest point(s). Note that the distance is converted in function from radians to kilometers - by multiplying by the radius of the earth and dividing the resulting meters by 1000. + by multiplying by the radius of the earth in km. """ X = coords_to_point_array([indexers[c] for c in self._index_coords]) @@ -201,7 +201,7 @@ def _query(self, indexers): concatenate=True, ) - return results, distances * 6371000 / 1000 + return results, distances * 6371 def _get_pos_indexers(self, indices, indexers): """Returns positional indexers based on the query results and the From e58322ee716c729bf4b424a8bffff6f3d35280c1 Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Mon, 23 Jan 2023 10:13:44 -0600 Subject: [PATCH 13/14] updates due to comments from benbovy --- src/xoak/accessor.py | 9 ++------- src/xoak/tests/test_accessor.py | 3 ++- 2 files changed, 4 insertions(+), 8 deletions(-) diff --git a/src/xoak/accessor.py b/src/xoak/accessor.py index f7c1294..cacae87 100644 --- a/src/xoak/accessor.py +++ b/src/xoak/accessor.py @@ -134,11 +134,7 @@ def index(self) -> Union[None, Index, Iterable[Index]]: return [wrp.index for wrp in index_wrappers] def _query(self, indexers): - """Find the distance(s) and indices of nearest point(s). - - Note that the distance is converted in function from radians to kilometers - by multiplying by the radius of the earth in km. - """ + """Find the distance(s) and indices of nearest point(s).""" X = coords_to_point_array([indexers[c] for c in self._index_coords]) if isinstance(X, np.ndarray) and isinstance(self._index, XoakIndexWrapper): @@ -201,7 +197,7 @@ def _query(self, indexers): concatenate=True, ) - return results, distances * 6371 + return results, distances def _get_pos_indexers(self, indices, indexers): """Returns positional indexers based on the query results and the @@ -292,7 +288,6 @@ def sel( result = result.to_dataset() # use same dimensions as indexers attrs = { - 'units': 'km', 'long_name': 'Distance from location to nearest comparison point.', } diff --git a/src/xoak/tests/test_accessor.py b/src/xoak/tests/test_accessor.py index 6492c7e..a9400be 100644 --- a/src/xoak/tests/test_accessor.py +++ b/src/xoak/tests/test_accessor.py @@ -100,4 +100,5 @@ def test_distances(): ) assert isinstance(output, xr.Dataset) - assert np.allclose(output['distances'], [0, 55.59746332]) + # this distance is in radians + assert np.allclose(output['distances'], [0, 0.00872665]) From 40a9ba204f2f249f4d88e647bed896eecc509f27 Mon Sep 17 00:00:00 2001 From: Kristen Thyng Date: Mon, 23 Jan 2023 13:22:41 -0600 Subject: [PATCH 14/14] updates for willirath --- src/xoak/accessor.py | 4 +++- src/xoak/tests/conftest.py | 1 + src/xoak/tests/test_accessor.py | 2 +- 3 files changed, 5 insertions(+), 2 deletions(-) diff --git a/src/xoak/accessor.py b/src/xoak/accessor.py index cacae87..919f757 100644 --- a/src/xoak/accessor.py +++ b/src/xoak/accessor.py @@ -281,7 +281,9 @@ def sel( result = self._xarray_obj.isel(indexers=pos_indexers) # save distances as a new variable in xarray object if name is input - if distances_name is not None: + if distances_name is not None or ( + isinstance(distances_name, str) and len(distances_name) == 0 + ): # need to have a Dataset instead of DataArray to add a new variable # otherwise goes in as a coordinate if not isinstance(result, xr.Dataset): diff --git a/src/xoak/tests/conftest.py b/src/xoak/tests/conftest.py index 3774cba..eec5782 100644 --- a/src/xoak/tests/conftest.py +++ b/src/xoak/tests/conftest.py @@ -1,3 +1,4 @@ +import dask import dask.array import numpy as np import pytest diff --git a/src/xoak/tests/test_accessor.py b/src/xoak/tests/test_accessor.py index a9400be..a35aae3 100644 --- a/src/xoak/tests/test_accessor.py +++ b/src/xoak/tests/test_accessor.py @@ -101,4 +101,4 @@ def test_distances(): assert isinstance(output, xr.Dataset) # this distance is in radians - assert np.allclose(output['distances'], [0, 0.00872665]) + np.testing.assert_allclose(output['distances'], [0, 0.008726646259971648])