diff --git a/notebooks/working_notebook.ipynb b/notebooks/working_notebook.ipynb new file mode 100644 index 0000000..e40434e --- /dev/null +++ b/notebooks/working_notebook.ipynb @@ -0,0 +1,1397 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### First implementation of RGDR's dbscan-based clustering & dimensionality reduction\n", + "This notebook outlines the current status of the RGDR implementation." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First we will load in some example data, and resample them using the `AdventCalendar`" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import xarray as xr\n", + "import s2spy.time\n", + "import s2spy.rgdr\n", + "\n", + "file_path = '../tests/test_rgdr/test_data'\n", + "field = xr.open_dataset(f'{file_path}/sst_daily_1979-2018_5deg_Pacific_175_240E_25_50N.nc')\n", + "target = xr.open_dataset(f'{file_path}/tf5_nc5_dendo_80d77.nc')\n", + "\n", + "cal = s2spy.time.AdventCalendar((8, 31), freq = \"30d\")\n", + "field_resampled = cal.resample(field)\n", + "target_resampled = cal.resample(target)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The target timeseries comes from already pre-clustered land surface temperature data." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "target_timeseries = target_resampled.sel(cluster=3).ts" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The precursor field consists of sea surface temperature data, with latitude and longitude dimensions:" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "With `_rgdr.correlation` we can determine correlation coefficient and p-values" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "field_resampled['corr'], field_resampled['p_val'] = (\n", + " s2spy.rgdr._rgdr.correlation(field_resampled.sst,\n", + " target_timeseries.sel(i_interval=0),\n", + " corr_dim='anchor_year')\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:      (anchor_year: 39, i_interval: 12, latitude: 5, longitude: 13)\n",
+       "Coordinates:\n",
+       "    index        (anchor_year, i_interval) int64 0 1 2 3 4 ... 464 465 466 467\n",
+       "    interval     (anchor_year, i_interval) object (1980-08-01, 1980-08-31] .....\n",
+       "  * latitude     (latitude) float64 47.5 42.5 37.5 32.5 27.5\n",
+       "  * longitude    (longitude) float64 177.5 182.5 187.5 ... 227.5 232.5 237.5\n",
+       "  * anchor_year  (anchor_year) int64 1980 1981 1982 1983 ... 2015 2016 2017 2018\n",
+       "  * i_interval   (i_interval) int64 0 1 2 3 4 5 6 7 8 9 10 11\n",
+       "    target       (i_interval) bool True False False False ... False False False\n",
+       "    tfreq        int64 5\n",
+       "    n_clusters   int64 6\n",
+       "    cluster      int64 3\n",
+       "Data variables:\n",
+       "    sst          (anchor_year, i_interval, latitude, longitude) float64 284.2...\n",
+       "    corr         (i_interval, latitude, longitude) float64 -0.08036 ... -0.00...\n",
+       "    p_val        (i_interval, latitude, longitude) float64 0.6267 ... 0.9864
" + ], + "text/plain": [ + "\n", + "Dimensions: (anchor_year: 39, i_interval: 12, latitude: 5, longitude: 13)\n", + "Coordinates:\n", + " index (anchor_year, i_interval) int64 0 1 2 3 4 ... 464 465 466 467\n", + " interval (anchor_year, i_interval) object (1980-08-01, 1980-08-31] .....\n", + " * latitude (latitude) float64 47.5 42.5 37.5 32.5 27.5\n", + " * longitude (longitude) float64 177.5 182.5 187.5 ... 227.5 232.5 237.5\n", + " * anchor_year (anchor_year) int64 1980 1981 1982 1983 ... 2015 2016 2017 2018\n", + " * i_interval (i_interval) int64 0 1 2 3 4 5 6 7 8 9 10 11\n", + " target (i_interval) bool True False False False ... False False False\n", + " tfreq int64 5\n", + " n_clusters int64 6\n", + " cluster int64 3\n", + "Data variables:\n", + " sst (anchor_year, i_interval, latitude, longitude) float64 284.2...\n", + " corr (i_interval, latitude, longitude) float64 -0.08036 ... -0.00...\n", + " p_val (i_interval, latitude, longitude) float64 0.6267 ... 0.9864" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "field_resampled" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Plot `lag 1` correlation map:" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhAAAADgCAYAAABSDV+5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkH0lEQVR4nO3debgcVZ3/8fcnNwlrMMRADEkEFBQBMWAEZ0REBAQXYNwAtzDiRGX8DS6jgMwwCuKAG4yjo0ZUoqKIKBoREAQiLmwJhLCEJUo0CSExbCaIkNz7/f1R55JO0+vtrq57uz+v56knVaeq+pxT3bn97XNOnVJEYGZmZtaMUUUXwMzMzEYeBxBmZmbWNAcQZmZm1jQHEGZmZtY0BxBmZmbWNAcQZmZm1jQHEMOcpDslHVjnmE9IOq8zJWqNpHmS3lt0Oayyot4fSR+QtErSOknP7nT+I42kAyUtL7oc1tscQAxzEbFHRMyrc8xnIqKhP/qSPinpe20p3DAh6dWSrpX0mKSlTZ57nKTf5lS0QvOXtFTSE+lLeXDZIY+8mixXSNqlZHsM8EXg0IjYOiIeKq501aX3qr/seh5YdLnMiuIAwpoiaXTRZajgceBbwMc6nfEwvR6l3pi+lAeXB4ouUAWTgM2BOyvtHGbX+Pqy6zmv6AKZFcUBxDCXfkUeXOeYp1sVJO2UfuHNlPRnSWsknZr2HQZ8Ajg6/Xq6LaU/S9I3Ja2UtELSpyX1pX3HSfqdpHMkPQScIelRSXuW5L9d+qW7vaRtJV0q6S+SHknrU3O6PABExE0R8V3gj82cJ+lFwNeAf0jX49GU/npJt0r6q6Rlkj5Zcs7g9T1e0p+BayT1SfpCutb3S/pgOmZ0Oqfi9a2Wf56aeX8k7SLp16llZ42kH5bs203SVZIelnSPpLc1mP91afW2VOeTgHtS2qOSrknHhaR/lXQfcF9Ke4Okhenz93tJe5W87t6SbpG0VtIPJV0o6dPNX6F8DLY0Sfp8uu73Szq8gfMmSPq2pAfSeT/tQHHNGuIAonvtD7wQeA1wmqQXRcQVwGeAH6ZfTy9Jx54PbAB2AfYGDgVKu0T2I/tyngScDvwEOLZk/9uAX0fEarLP1LeBHYHnAk8AX26kwJLenr4cqi3Pbf4yVBcRi4H3s/FX5fi063Hg3cB44PXAByQdVXb6q4AXAa8F/gU4HJgO7AOUH3s+Fa5vjfw3Ien/alyTRU1Wu5n35wzgSmBbYCrwv6k8WwFXAd8HtgeOAf5P0u71Mo+IA9LqS1Kdzwb2SGnjI+KgksOPIvvs7S5pb7JWpvcBzwa+DsyVtJmkscBPge8CE4AfAW+uVgZJ+9f5nO1fowp7p2DqXkn/qeZaR/YjC5YmAp8FvilJdc75LrAl2TXaHjinifzM8hURXobxAiwFDq5zzCeB76X1nYAAppbsvwk4pvzYtD0JeBLYoiTtWODatH4c8Oey/A4G/lCy/Tvg3VXKNh14pGR7HtmXZx7X6mBgaZPnHAf8ts4x5wLnlF3f55XsvwZ4X1k5Ahjd4PWtmX+Ln511wKNp+Wkz7w/wHWB26WcppR8N/KYs7evAfzVYrgB2KdkevKajy445qGT7q8AZZa9zD1kgdwDwAKCSfb8HPt3m6/k8YGeyIOzFwF3AKU18zpaUbG+Z6vicGudMBgaAbSvsOxBYnsfnxouXRpfh1Ldo7fVgyfrfgK2rHLcjMAZYWfJjaBSwrOSYZWXnXAtsKWk/YBXZl9AlAJK2JPuVdBjZL1eAcZL6IqJ/SDXpsFSvs4A9gbHAZmS/akuVXpMdqH69Grm+eToqIn41uNHk+/NxslaImyQ9AnwhIr5FVqf9yrpcRpP9Wm6n8us4U9L/K0kbS3btA1gREaVPBvxTm8tCRJR2kd0u6XSycTf/3eBLPP1/MiL+lj4P1f5fAkwDHo6IR5otq1knOIDoPeWPX11G9gt5YkRsaOSciOiXdBHZL+lVwKURsTbt/ihZ18l+EfGgpOnArUC9plokvYPsl2w1u0fEn+u9TpMqPY72+2TN+odHxN8lnUvW7FztvJVkTfyDppWs17u+dR+HK+lrwDur7P5TROxRZV8lDb8/EfEgWfcMqVn/V2kMwzKyLqtDmsh3KEqvzTLgzIg4s/wgSa8CpkhSSRDxXOAPlV5U0iuBy2vke3hE/KbB8tX9XLdgGTBB0viIeDTHfMyGxGMges8qYCdJowAiYiVZP/cXJG0jaZSk56c/yrV8n6wp+x1pfdA4sn71RyVNAP6r0YJFxAWx6Qj38qVi8JDKvDnZL31J2jz1iw/un6eSgZBlVgFTS49PdXg4BQ/7Am+vU/SLgBMlTZE0HjippE71rm+l/Muvy/trXJNmgofBujX0/kh6qzYOsHyE7AtzALgUeIGkd0kak5aXKRsUOjhgcGmNMqwi6w5oxjeA90vaT5mtlA12HQdcTzbG5N9SWd4E7FvthSLiN3U+ZxWDB0mHS5qU1ncD/hP4Wcn+Wp+zpqXPzuVk40u2TXU7oN55Zp3iAKL3DDbFPyTplrT+brLm4LvIviguJut/rSoibiQbbLgDm/6aOxfYAlgD3ABc0a6C13AA2ZfiZWwcGHhlyf5pZOM0KrmG7PbBByWtSWknAKdLWgucRhYg1PKNlN8isl/zl5F9oQ12CdS6vpXyz9O5NP7+vAy4UdI6YC5wYkT8MbU2HUo2ePIBsqb5s8m6eqD29YZsHM6cNGCxobs3ImI+WWvIl8mu4RKycQVExFPAm9L2w2SB7U8aed0mvQZYJOlxsvf4J2SDkgfVq/dQvAtYD9wNrAY+VOkgSZdL+kTJ9rrU0oKkV6b30KyttGm3oVl3Sb+gL4qIf+xgnocDX4uIHTuV53Ai6UqyYGNxgWU4n2yQ4X90KL+Of87MiuYxENbVImI5kOsfdUlbAK8ma4WYRNYtcEmeeQ5nEXFo0WXotE58zsyGm1y7MJRNgnS7sslf5qe0CcomoLkv/bttvdexp5so11VYPlH/bMuZgE+RNa3fCiwm6/owq6vK/+unuyDMhqtcuzDSQKoZEbGmJO2zZAPUzpJ0Mtk9zidVew0zMzMbfooYRHkkMCetz+GZs/aZmZnZMJd3ABHAlZIWSJqV0ial25MgG709KecymJmZWZvlPYhy/4hYIWl74CpJd5fujIiQVLEPJQUcswA0duxLx0zaPueiVpHnNDF1jPp7cXmP/ntxd+eMWl/UhJUFvtkFZs1T1eYP63KjirvoMba48evrtyru7v09phb3e3HBggVrImK7PPN47au3jDUPD1Tdf8uiJ38ZEYflWYZOyvVTHBEr0r+rJV1CNrnLKkmTI2KlpMlk9zZXOnc22Tz8bPbcaTH5pA/lWdTq+or7Ih13b19heU+4e31heW+x/K/FZDyquD+sMaa493rUslWF5U3dZ0nlaLOqc3flbv2OuX6P1bRqxhaF5T3/Cx8uLG9JbZ/evNyah/v5/RVTqu7ffIf7y2e0HdFy+4uZZoobN7hONvHMHWQT0sxMh82kZCY3MzOzkSqADfRXXbpNni0Qk4BL0gNjRgPfj4grJN0MXCTpeLIH3jQ0E52ZmdlwFgT9PTQ5Y24BRHpy3UsqpD9ENiWsmZlZ1whgPdXHQDRC0mHA/wB9wHkRcVaV495MNi3+y9JU7x3nmSjNzMzaIID1MfQAQlIf8BXgEGA5cLOkuRFxV9lx44ATgRuHXtrW+WFaZmZmbTJQY2nAvsCS9NC6p4ALyeZOKncG2QPsCrxXzwGEmZlZW0QET9VYgImS5pcss8peYgqwrGR7eUp7mqR9gGkR8YtcK9MAd2GYmZm1QVC3pWFNRMwY6utLGgV8kfQo+6I5gDAzM2uDQKyPluY2WQFMK9memtIGjQP2BOalOxyfA8yVdEQRAykdQJiZmbVJf2vTyt4M7CppZ7LA4Rjg7YM7I+Ix4OnJqCTNA/7dd2GYmZmNYNldGEMfWhgRGyR9EPgl2W2c34qIOyWdDsyPiLntKWl7OIAwMzNrgwHEU7Q2LX1EXAZcVpZ2WpVjD2wpsxY5gDAzM2uTgdbGQIwoDiDMzMzaIBBPRXEPxus0BxBmZmZtkN3G2TvTKzmAMDMza4MIt0CYmZnZEAy0dhvniOIAwszMrA2yMRC987XaOzU1MzPLUTYPhLswzMzMrAmB6PcgSjMzM2tG1gLRO1+rvVNTMzOzHAWi3xNJmZmZWTMi3AJhZmZmTZNv4zQzM7PmBPg2TjMzM2tOIN/GaWZmZs0JYCB8G6eZmZk1oddaIHIPlST1SbpV0qVp+3xJ90tamJbpeZfBzMysE/pR1aXbdKIF4kRgMbBNSdrHIuLiDuRtZmbWERFi/UDvNOzn2gIhaSrweuC8PPMxMzMrWpA9jbPa0m3yDpXOBT4OjCtLP1PSacDVwMkR8WT5iZJmAbMARo/fltHrihmY0r/lQCH5Aqzbubi81+5S3ECgvicmFJLvwGZRSL4AWzxQ3PWeeMeWheW91R2rCst7w9I/FZb3mMefKCzv/le8sLC8u10g1g94DETLJL0BWB0RC8p2nQLsBrwMmACcVOn8iJgdETMiYkbfVlvlVUwzM7O2GBxEWW3pNnn+7HkFcISkpcCFwEGSvhcRKyPzJPBtYN8cy2BmZtYxA4yqunSb3GoUEadExNSI2Ak4BrgmIt4paTKAJAFHAXfkVQYzM7NOiYD1A6OqLt2miOGiF0jaDhCwEHh/AWUwMzNrq0CeSKrdImIeMC+tH9SJPM3MzDopgPUOIMzMzKw5boEwMzOzJkW4BcLMzMyaFIgNngfCzMzMmtXqTJSSDpN0j6Qlkk6usP8jku6StEjS1ZJ2bHslGuQAwszMrA0C2DDQV3WpR1If8BXgcGB34FhJu5cddiswIyL2Ai4GPtveWjTOAYSZmVk7hBiosTRgX2BJRPwxIp4im4TxyE2yiLg2Iv6WNm8Apra1Dk3wGAgzM7M2CGBD7UGUEyXNL9meHRGzS7anAMtKtpcD+9V4veOBy5stZ7s4gDAzM2uDgHotDWsiYkY78pL0TmAG8Kp2vN5QOIAwMzNrg+wujJZGBqwAppVsT01pm5B0MHAq8KpKT7PuFAcQZmZm7RB1uzDquRnYVdLOZIHDMcDbSw+QtDfwdeCwiFjdSmatcgBhZmbWBg10YdQ+P2KDpA8CvwT6gG9FxJ2STgfmR8Rc4HPA1sCPsmdS8ueIOKLlwg+BAwgzM7M2aEMXBhFxGXBZWdppJesHt5RBGzmAMDMza5NooQVipHEAYWZm1gbR+hiIEcUBhJmZWZu4BcLMzMyaJPpbHAMxkjiAMDMza4NW78IYaRxAmJmZtUNAvwMIMzMza0a4C8PMzMyGIqLoEnSOAwgzM7M2iIABt0CYmZlZszyI0szMzJo2MNA7AUTubS2S+iTdKunStL2zpBslLZH0Q0lj8y6DmZlZ3gIRUX3pNp3orDkRWFyyfTZwTkTsAjwCHN+BMpiZmeUrsi6Maku3yTWAkDQVeD1wXtoWcBBwcTpkDnBUnmUwMzPrlBhQ1aXb5D0G4lzg48C4tP1s4NGI2JC2lwNTci6DmZlZR/g2zjKSXgB8FZgUEXtK2gs4IiI+XeOcNwCrI2KBpAObLZikWcAsgL4J49mwTX+zL9EW46asLSRfgLWrti4s73GT1hWW9wsm/qWQfBfcu2Mh+QI8OWFMYXmPWbuh/kF5eeKJ4vIu0qjifo0+a+lAYXl3uwiIHrqNs9GafgM4BVgPEBGLgGPqnPMK4AhJS4ELybou/gcYL2kwcJkKrKh0ckTMjogZETGjb+vivkjNzMwaFVF96TaNBhBbRsRNZWk1f7ZExCkRMTUidiILNq6JiHcA1wJvSYfNBH7WRHnNzMyGqerjH7pxDESjAcQaSc8ne9gYkt4CrBxinicBH5G0hGxMxDeH+DpmZmbDS9RYukyjgyj/FZgN7CZpBXA/8M5GM4mIecC8tP5HYN+mSmlmZjbcBV3Z0lBNQwFE+tI/WNJWwKiIKG5koZmZ2XDVhfM9VFMzgJD0kSrpAETEF3Mok5mZ2cjUhV0V1dRrgRicv+GFwMuAuWn7jUD5oEozM7Pe5S6MjSLiUwCSrgP2Gey6kPRJ4Be5l87MzGwk6aEWiEbvwpgEPFWy/VRKMzMzs0QDqroMN+lhl58f6vmN3oXxHeAmSZek7aPInmNhZmZmMOJu14yIfkn7D/X8Ru/COFPS5cArU9I/R8StQ83UzMys+wiGYUtDHbdKmgv8CHh8MDEiflLvxEafhfFcYA1wSWlaRPy5+bKamZl1qZH3qJHNgYfIHjcxKID2BBBkAyYHG2a2AHYG7gH2aLyMZmZmXSxoeR4ISYeRPTeqDzgvIs4q278Z2bCCl5J98R8dEUuHmFcf8FBE/PtQzm+0C+PFZZnuA5wwlAzNzMy6lVpogUhf6F8BDgGWAzdLmhsRd5UcdjzwSETsIukY4Gzg6KHkl8ZAvGKo5W20BaI801sk7TfUTM3MzOwZ9gWWpNmfkXQhcCRQGkAcCXwyrV8MfFmSIob8vM+FeY+BKJ2RchSwD/BAk4U0MzPranVu15woaX7J9uyImF2yPQVYVrK9HCj/sf70MRGxQdJjZA+mXDPEIuc+BmJcyfoGsjERP260dGZmZl2v/m2cayJiRmcK05iI+OehnttoAHFXRPyoNEHSW8maPMzMzIzWxkAAK4BpJdtTU1qlY5ZLGg08i6wFYUgkTQX+FxgcC/Eb4MSIWF7v3EZnojylwTQzM7PeNVBjqe9mYFdJO0saCxzDxmdQDZoLzEzrbwGuaWH8A8C302vukJafp7S66j2N83DgdcAUSV8q2bUNWVeGmZmZAYpsGao0puGDwC/JbuP8VkTcKel0YH5EzAW+CXxX0hLgYbIgoxXbRURpwHC+pA81cmK9LowHgPnAEcCCkvS1wIebKaGZmVnXa3Emyoi4DLisLO20kvW/A29tKZNNPSTpncAP0vaxNNglUu9pnLcBt0m6ICLc4mBmZlZDKy0QBXkP2RiIc8iGgP4eOK6RE+t1YVwUEW8jmyv7GZclIvZquqhmZmbdKFoeRFmE04GZEfEIgKQJwOfJAoua6nVhnJj+fUNLxTMzM+sFI68FYq/B4AEgIh6WtHcjJ9a8CyMiVqbVEyLiT6ULnsrazMxsExqovgxToyRtO7iRWiAamuKh0ds4D6mQdniD55qZmfWGqLEMT18Arpd0hqQzyMZAfLaRE+uNgfgAWUvD8yQtKtk1DvjdEAtrZmbWfVq8jbMIEfGdNL324FTWbyp7eFdV9Zopvg9cDvw3cHJJ+tqIeLjpkpqZmXWz4dtVUVUKGBoKGkrVGwPxWEQsjYhj07iHJ8gaYraW9Nxa50raXNJNkm6TdKekT6X08yXdL2lhWqY3W2gzM7PhRmycTKrS0m0afRrnG4Evkk1zuRrYEVgM7FHjtCeBgyJinaQxwG8lXZ72fSwiLh56sc3MzIaZkXkb55A1Oojy08DLgXsjYmfgNcANtU6IzLq0OSYtXRiDmZmZJSNvEOWQNRpArI+Ih8hu9xgVEdcCdR9JKqlP0kKyVourIuLGtOtMSYsknSNpsyGV3MzMbJgZgbdxDlmjj/N+VNLWwHXABZJWA4/XOyki+oHpksYDl0jak+wpng8CY4HZwElkM2FtQtIsYBbA6InPYvSz/95gUdtrvx3+VEi+AG968fzC8u4rcCTQdWt3KyTfO7Z5TiH5ArBsbGFZr5tWXAy/+dY7FZb32IcnF5Y3qx4rLOst/rK+sLxfc+BnCsu7I4IROYhyqBptgTiSbADlh4ErgD8Ab2w0k4h4FLgWOCwiVqbujSfJHhm6b5VzZkfEjIiY0bfNVo1mZWZmVhgPoiwTEaWtDXMaOUfSdmRdH49K2oJsMqqzJU2OiJWSBBwF3NFkmc3MzIalbuyqqKbeRFJrqTz0Q2TjJLepcfpkYI6kPrKWjosi4lJJ16TgQsBC4P1DKrmZmdlw04UtDdXUe5z3uKG+cEQsAp7xQI6IOKjC4WZmZiNat3ZVVNPoIEozMzOrxwGEmZmZNctjIMzMzKw5PTYTpQMIMzOzdnEXhpmZmTXLLRBmZmbWNN+FYWZmZs3psamsHUCYmZm1gXALhJmZmQ2BBnongnAAYWZm1g6+jdPMzMyGpHcaIBxAmJmZtUsvtUCMKroAZmZmXSE2PlCr0tIKSRMkXSXpvvTvthWOmS7pekl3Slok6ejWcq3NAYSZmVkbiKwFotrSopOBqyNiV+DqtF3ub8C7I2IP4DDgXEnjW865CgcQZmZm7RJRfWnNkcCctD4HOOqZWce9EXFfWn8AWA1s12rG1XgMhJmZWTvUvwtjoqT5JduzI2J2g68+KSJWpvUHgUm1Dpa0LzAW+EODr980BxBmZmZtov6au9dExIyq50q/Ap5TYdeppRsREVL1URWSJgPfBWZGRG7DOh1AmJmZtUkrgyUj4uCqryutkjQ5IlamAGF1leO2AX4BnBoRNwy9NPV5DISZmVk7RDYTZbWlRXOBmWl9JvCz8gMkjQUuAb4TERe3mmE9DiDMzMzaJWosrTkLOETSfcDBaRtJMySdl455G3AAcJykhWmZ3nLOVbgLw8zMrA0UbWlpqCgiHgJeUyF9PvDetP494Hu5FKACBxBmZmZt4qdxmpmZWdN6aSprBxBmZmbtEEB/7zRB5DaIUtLmkm6SdFual/tTKX1nSTdKWiLph2nUqJmZ2YiX17MwhqM878J4EjgoIl4CTAcOk/Ry4GzgnIjYBXgEOD7HMpiZmXVMjrdxDju5BRCRWZc2x6QlgIOAwftTK87nbWZmNuLUuoWz++KHfMdASOoDFgC7AF8hm5P70YjYkA5ZDkypcu4sYBbAtCl93PbKb+RZ1Kp+/fdxheRbtC8tqzohWu7uuG9qIflus3hMIfkCbP5wcX9dNl+zof5BOdny7oqT6XXEwOq/FJY3E57xJOaOGXvTvYXlzZjuHnYnQB4D0R4R0R8R04GpwL7Abk2cOzsiZkTEjInP7suriGZmZm2jiKpLt+lIOBgRj0q6FvgHYLyk0akVYiqwohNlMDMzy1UEdOFYh2ryvAtjO0nj0/oWwCHAYuBa4C3psIrzeZuZmY1EvTSIMs8WiMnAnDQOYhRwUURcKuku4EJJnwZuBb6ZYxnMzMw6IzyRVFtExCJg7wrpfyQbD2FmZtZdurCloZruHhJrZmbWQd04WLIaBxBmZmbt0GNTWTuAMDMzawPRnbdrVuMAwszMrF0GemcUpQMIMzOzdgigd+IHBxBmZmbtIrdAmJmZWVMi3IVhZmZmQ9A78YMDCDMzs3ZxF4aZmZk1J/BMlGZmZtYsj4EwMzOzofBEUmZmZtaUCOjvL7oUHeMAwszMrB0C6HcXhpmZmTXLXRhmZmbWnN4aRDmq6AKYmZl1hSALIKotLZA0QdJVku5L/25b49htJC2X9OWWMq3DAYSZmVm75BRAACcDV0fErsDVabuaM4DrWs2wHgcQZmZmbRHZRFLVltYcCcxJ63OAoyodJOmlwCTgylYzrMdjIMzMzNohIGrfxjlR0vyS7dkRMbvBV58UESvT+oNkQcImJI0CvgC8Ezi4wdcdMgcQZmZm7VB/Hog1ETGj2k5JvwKeU2HXqZtmEyGpUpPGCcBlEbFcUiMlbokDCDMzs3Zp4TbOiKjaaiBplaTJEbFS0mRgdYXD/gF4paQTgK2BsZLWRUSt8RJD5gDCzMysLaJeF0Yr5gIzgbPSvz97Ru4R7xhcl3QcMCOv4AFyHEQpaZqkayXdJelOSSem9E9KWiFpYVpel1cZzMzMOmbwaZz5DKI8CzhE0n1k4xvOApA0Q9J5rb74UOTZArEB+GhE3CJpHLBA0lVp3zkR8fkc8zYzM+uooO4gyqG/dsRDwGsqpM8H3lsh/Xzg/FwKk+QWQKTRoivT+lpJi4EpeeVnZmZWqAgIz0TZVpJ2AvYGbkxJH5S0SNK3as2mZWZmNpJEf3/Vpdsocn7wh6StgV8DZ0bETyRNAtaQtfacAUyOiPdUOG8WMCttvhC4J9eCVjeRrLy9phfr3Yt1Bte71/RqvV8YEePyzEDSFWTXt5o1EXFYnmXopFwDCEljgEuBX0bEFyvs3wm4NCL2zK0QLZI0v9Z9u92qF+vdi3UG17vocnSa623tkuddGAK+CSwuDR7S/auD/gm4I68ymJmZWT7yvAvjFcC7gNslLUxpnwCOlTSdrAtjKfC+HMtgZmZmOcjzLozfApXm0rwsrzxz0ug85d2mF+vdi3UG17vXuN7WFrkPojQzM7Pu48d5m5mZWdN6PoBIc1GslnRHSdp0STekqbbnS9o3pUvSlyQtSfNY7FNcyYeuSp1fIul6SbdL+rmkbUr2nZLqfI+k1xZT6tbVmF59gqSrJN2X/t02pXfL+12t3m9N2wOSZpSdM6Lf8xp1/pyku9P7eYmk8SXnjOg6Q816n5HqvFDSlZJ2SOld/Rkv2f9RSSFpYtruinoXLiJ6egEOAPYB7ihJuxI4PK2/DphXsn452diOlwM3Fl3+Ntb5ZuBVaf09wBlpfXfgNmAzYGfgD0Bf0XUYYr0nA/uk9XHAval+nwVOTuknA2d32ftdrd4vIptjZR7ZQ3cGjx/x73mNOh8KjE7pZ5e81yO+znXqvU3JMf8GfC2td/VnPG1PA34J/AmY2E31Lnrp+RaIiLgOeLg8GRj8Bf4s4IG0fiTwncjcAIwvuy11RKhS5xcA16X1q4A3p/UjgQsj4smIuB9YAuzbkYK2WUSsjIhb0vpaYHB69SOBOemwOcBRab1b3u+K9Y6IxRFRaYK2Ef+e16jzlRGxIR12AzA1rY/4OkPNev+15LCtyP7GQZd/xtPuc4CPs7HO0CX1LlrPBxBVfAj4nKRlwOeBU1L6FGBZyXHL6Z7ne9xJ9p8K4K1kUTt0aZ216fTqkyJ7dgvAg8CktN51ddczp5WvpKvqXaPO7yH7FQpdVmd4Zr0lnZn+pr0DOC0d1tX1lnQksCIibis7rOvqXQQHEJV9APhwREwDPkw2IVa3ew9wgqQFZE2ATxVcntwom179x8CHyn6ZERHBpr9UukateneranWWdCrZE4MvKKpseapU74g4Nf1NuwD4YJHly0tpvcne30+wMViyNnMAUdlM4Cdp/UdsbMpcwcZf5pA1f67oYLlyExF3R8ShEfFS4AdkfcDQZXVWNr36j4ELImLwPV412HyZ/l2d0rum7lXqXU1X1LtanSUdB7wBeEcKGKFL6gwNvdcXsLGLspvr/Xyy8Sy3SVpKVrdbJD2HLqp3kRxAVPYA8Kq0fhBwX1qfC7w7jeB9OfBYSdP3iCZp+/TvKOA/gK+lXXOBYyRtJmlnYFfgpmJK2Rqp8vTqZHWcmdZnAj8rSR/x73eNelcz4t/zanWWdBhZf/gREfG3klNGfJ2hZr13LTnsSODutN61n/GIuD0ito+InSJiJ7Juin0i4kG6pN6FK3oUZ9EL2a/tlcB6sg/Y8cD+wAKyUdk3Ai9Nxwr4Ctmv89spGbk+kpYqdT6RbOTyvcBZpEnG0vGnpjrfQ7o7ZSQu6X0NYBGwMC2vA54NXE0WKP4KmNBl73e1ev9Tev+fBFaRPfSuK97zGnVeQtb3PZj2tW6pc516/5jsuUOLgJ+TDazs+s942TFL2XgXRlfUu+jFM1GamZlZ09yFYWZmZk1zAGFmZmZNcwBhZmZmTXMAYWZmZk1zAGFmZmZNcwBh1iGS1uXwmkdIOjmtHyVp9yG8xjyVPY3TzKweBxBmI1hEzI2Is9LmUWRPXjQzy50DCLMOS7PffU7SHZJul3R0Sj8wtQZcLOluSRekGfaQ9LqUtkDSlyRdmtKPk/RlSf8IHEH2ELiFkp5f2rIgaWKazhdJW0i6UNJiSZcAW5SU7VBJ10u6RdKP0rMFzMyeYXTRBTDrQW8CpgMvASYCN0safJT63sAeZNOp/w54haT5wNeBAyLifkk/KH/BiPi9pLnApRFxMUCKPSr5APC3iHiRpL2AW9LxE8mmMT84Ih6XdBLwEeD0NtTZzLqMAwizztsf+EFE9JM9yOvXwMuAvwI3RcRyAEkLgZ2AdcAfI+L+dP4PgFkt5H8A8CWAiFgkaVFKfzlZF8jvUvAxFri+hXzMrIs5gDAbXp4sWe+ntf+jG9jYTbl5A8cLuCoijm0hTzPrER4DYdZ5vwGOltQnaTuyFoFaT368B3iepJ3S9tFVjlsLjCvZXgq8NK2/pST9OuDtAJL2BPZK6TeQdZnskvZtJekFjVTIzHqPAwizzruE7KmBtwHXAB+P7BHDFUXEE8AJwBWSFpAFCo9VOPRC4GOSbpX0fODzwAck3Uo21mLQV4GtJS0mG9+wIOXzF+A44AepW+N6YLdWKmpm3ctP4zQbASRtHRHr0l0ZXwHui4hzii6XmfUut0CYjQz/kgZV3gk8i+yuDDOzwrgFwszMzJrmFggzMzNrmgMIMzMza5oDCDMzM2uaAwgzMzNrmgMIMzMza5oDCDMzM2va/wf4UgzmciH+MQAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "field_resampled.corr.sel(i_interval=1).plot(cmap='viridis', size=3, aspect=3)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Next we can cluster the remaining grid points with `DBSCAN`.\n", + "\n", + "Compare these results with the [prototype-notebook](https://github.com/AI4S2S/s2spy/blob/6da3233168bc7b083b084a14f7afe5618b8c82c9/notebooks/prototype_RGDR.ipynb)." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhEAAADgCAYAAAC9zzSHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAoQ0lEQVR4nO3deZxcVZ3+8c+TsAooYCAgQUHAJTpsRmBGVEREQAVEURhUcGNAHWGcUUD9iYMyA24wjI4YEQUFEVE0g4CggLiwhS3syioJ24QdVCDJ8/vjniZFU1t3VXV1Vz/v1+u++m6n7jm3KqlvnXsW2SYiIiJipKb0OwMRERExMSWIiIiIiFFJEBERERGjkiAiIiIiRiVBRERERIxKgoiIiIgYlQQR45yk6yRt0+KcT0s6bmxy1BlJF0j6UL/zEfX16/2RtL+keyU9Jun5Y339iUbSNpLm9zsfEQkixjnbr7B9QYtz/sN2W//xS/q8pB90JXPjhKQ3SDpf0sOSbh9h2n0k/a5HWevr9SXdLumv5Yt5aHlBL641wnxZ0oY128sCXwO2t72y7fv7l7vGynu1eNj93Kbf+YropwQRMSKSlul3Hup4HDge+ORYX3ic3o9abytfzEPLXf3OUB3TgRWA6+odHGf3+KJh9/OCfmcoop8SRIxz5dfkdi3Oebp2QdJ65Zfe3pL+LGmhpM+UYzsAnwbeXX5FXV32P0/SdyTdLWmBpC9KmlqO7SPp95KOknQ/8AVJD0l6Zc311yi/eNeUtJqkMyT9n6QHy/qMHt0eAGxfavv7wK0jSSfp5cCxwN+X+/FQ2f8WSVdKekTSnZI+X5Nm6P5+UNKfgfMkTZX01XKvb5P0sXLOMiVN3fvb6Pq9NJL3R9KGkn5TangWSvpRzbGXSTpX0gOSbpL0rjavf2FZvbqU+SDgprLvIUnnlfMs6aOS/gT8qex7q6SryufvD5I2rnndzSRdIelRST+SdIqkL478DvXGUI2TpK+U+36bpB3bSLe6pO9Kuquk+9kYZDeibQkiBtfWwEuBNwKfk/Ry22cD/wH8qPyK2qSc+z1gEbAhsBmwPVD7eGRLqi/o6cBhwE+BPWuOvwv4je37qD5T3wVeBLwQ+Cvw9XYyLOkfyxdEo+WFI78Njdm+AdiPpb8uVy2HHgfeB6wKvAXYX9Kuw5K/Hng58Gbgw8COwKbA5sDwc79Hnfvb5PrPIOl/mtyTeSMs9kjeny8A5wCrATOA/y75WQk4FzgZWBPYA/gfSTNbXdz268rqJqXMRwKvKPtWtb1tzem7Un32ZkrajKq26Z+A5wPfAuZIWl7ScsDPgO8DqwM/Bt7RKA+Stm7xOdu6SRE2KwHVHyX9P42slmRLqoBpGvAl4DuS1CLN94HnUN2jNYGjRnC9iN6znWUcL8DtwHYtzvk88IOyvh5gYEbN8UuBPYafW7anA08AK9bs2xM4v6zvA/x52PW2A26p2f498L4GedsUeLBm+wKqL9Be3KvtgNtHmGYf4HctzjkaOGrY/X1xzfHzgH8alg8Dy7R5f5tev8PPzmPAQ2X52UjeH+BEYHbtZ6nsfzfw22H7vgUc2ma+DGxYsz10T5cZds62NdvfBL4w7HVuogrmXgfcBajm2B+AL3b5fr4YWJ8qEPs74HrgkBF8zm6u2X5OKeNaTdKsDSwBVqtzbBtgfi8+N1myjGQZT88ao7vuqVn/C7Byg/NeBCwL3F3zo2gKcGfNOXcOS3M+8BxJWwL3Un0RnQ4g6TlUv5Z2oPoFC7CKpKm2F4+qJGOslOsI4JXAcsDyVL9ua9XekxfQ+H61c397aVfbvxraGOH78ymq2ohLJT0IfNX28VRl2nLY45dlqH41d9Pw+7i3pH+u2bcc1b03sMB27WyCd3Q5L9iufVx2jaTDqNrh/GebL/H0v0nbfymfh0b/LgHWBR6w/eBI8xoxVhJETD7Dp229k+qX8jTbi9pJY3uxpFOpflHfC5xh+9Fy+F+pHqNsafseSZsCVwKtqm2RtBfVL9pGZtr+c6vXGaF609ieTFXFv6Ptv0k6mqoKulG6u6mq+4esW7Pe6v62nEZX0rHAexocvsP2Kxocq6ft98f2PVSPaihV/L8qbRrupHp89aYRXHc0au/NncDhtg8ffpKk1wPrSFJNIPFC4JZ6LyrptcBZTa67o+3ftpm/lp/rDtwJrC5pVdsP9fA6EaOWNhGTz73AepKmANi+m+q591clPVfSFEkblP+YmzmZqlp7r7I+ZBWq5+wPSVodOLTdjNk+yc9s+T58qRtAlDyvQPWLX5JWKM/Jh45foJrGkcPcC8yoPb+U4YESQGwB/GOLrJ8KHCBpHUmrAgfVlKnV/a13/eH3Zb8m92QkAcRQ2dp6fyTtrqWNLh+k+tJcApwBvETSeyUtW5ZXq2ooOtSI8PYmebiX6tHASHwb2E/SlqqspKoB7CrARVRtTj5e8rIbsEWjF7L92xafs7oBhKQdJU0v6y8D/h/w85rjzT5nI1Y+O2dRtTdZrZTtda3SRYylBBGTz1C1/P2Srijr76OqGr6e6sviNKrnsQ3ZvoSqAeILeOavuqOBFYGFwMXA2d3KeBOvo/piPJOljQXPqTm+LlW7jXrOo+paeI+khWXfR4DDJD0KfI4qSGjm2+V686h+1Z9J9aU29Hig2f2td/1eOpr2359XA5dIegyYAxxg+9ZS67Q9VYPKu6iq6Y+keuwDze83VO1yTiiNGNvq1WF7LlWtyNep7uHNVO0MsP0ksFvZfoAquP1pO687Qm8E5kl6nOo9/ilVQ+Uhrco9Gu8FngJuBO4DDqx3kqSzJH26ZvuxUuOCpNeW9zCi6/TMx4gRg6X8kj7V9j+M4TV3BI61/aKxuuZ4IukcqoDjhj7m4XtUDQ8/O0bXG/PPWcR4kDYRMdBszwd6+h+7pBWBN1DVRkynekRwei+vOZ7Z3r7feRhrY/E5ixiPevo4Q9VASdeoGiBmbtm3uqpBav5U/q7W6nXi6erKx+osn26dOnpMwL9TVbNfCdxA9RgkoqUG/66ffhwRk4ekdVUN4X+9qnmTDqhzjiQdI+lmSfMkbd6PvD6dn14+ziiNq2bZXliz70tUjdaOkHQwVR/ogxq9RkRExGQgaW1gbdtXlEbDl1N1076+5pydgH8GdqIawOy/bG/ZlwzTn4aVuwAnlPUTePbofhEREZOO7bttX1HWH6Wq1Vxn2Gm7ACe6cjGwagk++qLXQYSBcyRdLmnfsm966boEVavu6T3OQ0RExIQiaT2qYfIvGXZoHZ45ENt8nh1ojJleN6zc2vYCSWsC50q6sfagbUuq+zylBB37Ami55V617PQ1e5zViMlp+Tsf73cWJp0n1l2p31mYdJ68c/5C22v08hpvfsNzvPCBJQ2PXzHvieuAv9Xsmm179vDzJK0M/AQ40PYjXc9oF/U0iLC9oPy9T9LpVAPA3Ctpbdt3lyqY+xqknU01bj/Lv3Bdr/NvB/YyqxGT1gYHXtzvLEw6t/zbVv3OwqRz2wH/1vWh0Idb+MBi/nB240qBFV5w299sz2r2GpKWpQogTrJdb7yTBTxzVNwZZV9f9OxxRhlRbpWhdarBaa6lGrRm73La3tSM+BYRETFRGVjE4oZLK5IEfAe4wfbXGpw2B3hf6aWxFfBwTROBMdfLmojpwOnVPWEZ4GTbZ0u6DDhV0gepJslpa8S6iIiI8cyYxZ31eHwN1Sil10i6quz7NNVIvNg+lmq01J2oRm39C/D+Ti7YqZ4FEWXGu03q7L+favjYiIiIgWHgKRq3iWiZ3v4dLSZ1K5PMfXTUF+myjFgZERHRBQae8uiDiIkoQURERESXTK4QIkFEREREV9jmyUk2qWWCiIiIiC4wqYmIiIiIUTDiKTdtFzlwEkRERER0yeLmnSsGToKIiIiILqh6Z/RjXsv+SRARERHRBUsQTzK139kYUwkiIiIiumRJ2kRERETESBnxpFMTERERESNUdfFMm4iIiIgYIXvy1URMrpApIiKih5aghksrko6XdJ+kaxsc30bSw5KuKsvnul6AEUpNRERERBdUbSI6+lr9HvB14MQm5/zW9ls7uUg3JYiIiIjogmqciNE/zrB9oaT1upahMZDHGREREV1gxGKmNFyAaZLm1iz7juIyfy/paklnSXpFl4swYqmJiIiI6IKqJqLp1+pC27M6uMQVwItsPyZpJ+BnwEYdvF7HUhMRERHRBUYsduOl49e3H7H9WFk/E1hW0rSOX7gDqYmIiIjoArtlTURHJK0F3Gvbkragqgi4v2cXbEOCiIiIiK5orytnw9TSD4FtqNpOzAcOBZYFsH0s8E5gf0mLgL8Ce9h2p7nuRIKIiIiILjB01MXT9p4tjn+dqgvouJEgIiIioguMOuriOREliIiIiOgCA0s8uforJIiIiIjogslYE9HzkEnSVElXSjqjbH9P0m01Y39v2us8REREjIXFqOEyiMaiJuIA4AbguTX7Pmn7tDG4dkRExJiwxVNLJlcFf09rIiTNAN4CHNfL60RERPSb6WwWz4mo1yHT0cCngFWG7T+8TGH6a+Bg208MT1jGFN8XYOpqq/U4mxGT1y1Hb9W3a29w4MV9u3Y/9bPc/Xy/B50RTy1Jm4iukPRW4D7blw87dAjwMuDVwOrAQfXS255te5btWVNXXqlX2YyIiOiKoYaVjZZB1MvHGa8BdpZ0O3AKsK2kH9i+25UngO8CW/QwDxEREWNmCVMaLoOoZ6WyfYjtGbbXA/YAzrP9HklrA0gSsCtwba/yEBERMVZseGrJlIbLIOpHM9KTJK0BCLgK2K8PeYiIiOgqoww21Qu2LwAuKOvbjsU1IyIixpKBpzoIIiQdDwy1J3xlneMC/gvYCfgLsI/tK0Z9wS6YXCFTREREz1Q1EY2WNnwP2KHJ8R2BjcqyL/DNjrPcock1KkZERESP2J3VRNi+UNJ6TU7ZBTixTP99saRVJa1t++5RX7RDCSIiIiK6wIhFzceJmCZpbs32bNuzR3CJdYA7a7bnl30JIiIiIia6FiNTLrQ9a6zyMhYSRERERHSBoVVNRKcWAOvWbM8o+/omDSsjIiK6wWJJk6UL5gDvU2Ur4OF+toeA1ERERER0hYFFnXXx/CGwDVXbifnAocCyALaPBc6k6t55M1UXz/d3luPOJYiIiIjoAkNHNQ6292xx3MBHR32BHkgQERER0QVV74zJ1UogQUREREQ3uLPHGRNRgoiIiIgu6PRxxkSUICIiIqILBuFxhqQpwMq2H2nn/Ild2oiIiHHEVsNlvJJ0sqTnSloJuBa4XtIn20mbICIiIqILXNpENFrGsZml5mFX4CxgfeC97SQc16WKiIiYSCZiTQSwrKRlqYKIObafomri0VLaRERERHSFWDwx20R8C7gduBq4UNKLgLbaRCSIiIiI6IKJ2jvD9jHAMTW77pD0hnbSJoiIiIjoBsPiCRRESPpEi1O+1uo1EkRERER0gSfe44xVOn2BBBERERFd4raaIzYmaQfgv4CpwHG2jxh2fB/gyyydAvzrto8bzbVs/3sHWQXSOyMiIqIrbFiyZErDpRVJU4FvADsCM4E9Jc2sc+qPbG9allEFEMOu+xJJv5Z0bdneWNJn20mbICIiIqJLllgNlzZsAdxs+1bbTwKnALv0NMOVbwOHAE8B2J4H7NFOwgQRERERXbJkiRoubVgHuLNme37ZN9w7JM2TdJqkdbuQ7efYvnTYvkXtJOx5ECFpqqQrJZ1RtteXdImkmyX9SNJyvc5DREREr5nGA02VwaamSZpbs+w7isv8L7Ce7Y2Bc4ETupD1hZI2oAwwJemdwN3tJByLhpUHADcAzy3bRwJH2T5F0rHAB4FvjkE+IiIiesctx4lYaHtWk+MLgNqahRksbUBZXcK+v2bzOOBLI81mHR8FZgMvk7QAuA3Yq52EPa2JkDQDeAtVQZEkYFvgtHLKCVTDbEZEREx4XqKGSxsuAzYqNfbLUbVLmFN7gqS1azZ3pvqR3lmeqzYY2wFrAC+zvbXtO9pJ2+uaiKOBT7G0L+rzgYdsDz1rafS8JyIiYsLppIun7UWSPgb8kqqL5/G2r5N0GDDX9hzg45J2pmqz8ACwT6d5lvR84FBga8CSfgccNqzWo662gghJL6F65DDd9islbQzsbPuLTdK8FbjP9uWStmnnOsPS7wvsCzB1tdVGmjwiImJM2eAOB5uyfSZw5rB9n6tZP4SqJ0U3nQJcCLyjbO8F/AjYrlXCdks7mu4frwF2lnR7yeC2VANorCppKHh51vOeIbZn255le9bUlVdqM5sRERH9YzdexrG1bX/B9m1l+SIwvZ2E7QYRI+7+YfsQ2zNsr0cVcJxney/gfOCd5bS9gZ+3mYeIiIhxrHF7iDbbRPTLOZL2kDSlLO+ieqTSUrtBxKi7f9RxEPAJSTdTtZH4zihfJyIiYnxxk2WckfSopEeADwMnA0+W5RRKc4JW2m1YWa/7x3vazajtC4ALyvqtVKNyRUREDA4z3mscnsH22EzAVb74t5O0EjDF9qOdXjgiImLgTKCpwGtJWg3YCFhhaJ/tC1ulaxpENJprvBruAWy3nGs8IiJi0hiHjy1akfQhqoEhZwBXAVsBF1F1iGiqVZuIVcoyC9ifakyHdYD9gM1HneOIiIhB444Hm+qXA4BXA3fYfgOwGfBQOwmb1kQMzTUu6UJg86HHGJI+D/xi9PmNiIgYQBOwJgL4m+2/SULS8rZvlPTSdhK227ByOlWLzSFP0mYf0oiIiMlC47vGoZH5klYFfgacK+lBoKvDXp8IXCrp9LK9K92ZOSwiImIwjNOunK3YfntZ/byk84HnAWe3k7bd3hmHSzoLeG3Z9X7bV444pxEREQNLMIFqIiStXmf3NeXvylRzczTV7twZLwQWAqfX7rP953bSR0RETApL+p2BEbmcqu6kNvIZ2jbw4lYv0O7jjF+wtJJmRWB94CbgFe3mNCIiYqCZjseJkLQD1TxTU4HjbB8x7PjyVE0MXgXcD7zb9u2jyq69fpt5eoXt6+oda2vYa9t/Z3vjsmxENeLkRe1nNSIiYvBpSeOlZVppKvANYEdgJrCnpJnDTvsg8KDtDYGjgCO7W4K6vt/owKjmLLV9BbDlqLMTERERw20B3Gz7VttDc1jsMuycXVjaseE04I0aGgGydxq+frttImpHrpxCNdDUXR1mKiIiYqC06OI5TdLcmu3ZtmfXbK8D3FmzPZ9n/2B/+hzbiyQ9TDWZ5cJRZ7q1hn1O2m0TUTtJxyKqNhI/6SRHERERA6V1F8+FtmeNTWbGRrtBxPW2f1y7Q9LuwI8bnB8RETHptNP2oYkFwLo12zPKvnrnzJe0DNWYDveP9oLlUcgM23c2Oe3JRgfabRNxSJv7IiIiJq8lTZbWLgM2krS+pOWAPYA5w86ZA+xd1t8JnGd71ENclbRntjhnq0bHWs3iuSOwE7COpGNqDj2X6rFGREREAHK1jFZp4/Ax4JdUXTyPt32dpMOAubbnAN8Bvi/pZqrBoPboPOdcIenVti8bacJWjzPuAuYCO1MNSjHkUeBfRnqxiIiIgdbhiJW2z2RYzYDtz9Ws/w3YvaOLPNuWwF6S7gAepww2ZXvjVglbzeJ5NXC1pJNsp+YhIiKiiU5qIvrozaNN2Opxxqm23wVcKT371rQTpUREREwK7rhhZV/YvkPS1sBGtr8raQ2quTNaavU444Dy962dZDAiImJSmIA1EZIOBWYBLwW+CywL/AB4Tau0TXtn2L67rH7E9h21C/CRzrIdERExWDoZ9rqP3k7V9vFxANt38czxoRpqt4vnm+rs27HNtBEREZODmyzj15Olq6cBJK3UbsJWbSL2p6pxeLGkeTWHVgF+P4qMRkREDKYOu3j20amSvgWsKunDwAeA49pJ2KpNxMnAWcB/AgfX7H/U9gOjyWlERMTAGt+PLeqy/RVJbwIeoWoX8Tnb57aTtlUXz4eBh4E9ASStCawArCxpZdt/bpRW0grAhcDy5Tqn2T5U0veA15fXBdjH9lXtZDYiImK8EhOzJkLSkbYPAs6ts6+pttpESHqbpD8BtwG/AW6nqqFo5glgW9ubAJsCO0gaGjrzk7Y3LctV7eQhIiJiXPOEbVg56naP7Tas/CKwFfBH2+sDbwQubpbAlcfK5rJlmYAxWkRERJsmUMNKSftLugZ4qaR5NcttwLxW6aH9IOIp2/cDUyRNsX0+VZ/SVhmcKukq4D7gXNuXlEOHl4weJWn5NvMQERExrk2wmoiTgbdRTer1tprlVbbf084LtDsV+EOSVqZq43CSpPso/Umbsb0Y2FTSqsDpkl5JNfvnPcBywGzgIOCw4Wkl7QvsCzB1tdXazOZgee1W1/c7C33x24tn9jsLMUZuObrh5IADbYMDm1bkxkRletawUtLqwI+A9aiaFLzL9oN1zlsMXFM2/2x750avOdTuUdJngXtsPyFpG2BjSSfafqhVvtqtidgF+CvVpFtnA7dQRSttKRk5H9jB9t3lUccTVCNjbdEgzWzbs2zPmrpy211WIyIi+mZoJs96S4cOBn5teyPg1zyzx2Stv9a0OWwYQAzzE2CxpA2pftyvS1VL0VJbQYTtx20vtr3I9gm2jymPNxqStEapgUDSilQNN26UtHbZJ2BX4Np28hARETHe9fBxxi7ACWX9BKrvz25ZUibZ3A34b9ufBNZuJ2GrwaYepX5zkKFpQp/bJPnawAmSplIFK6faPkPSeWVyDwFXAfu1k9GIiIhxr3cNKKfXTEVxDzC9wXkrSJoLLAKOsP2zNl77KUl7Au9j6VOGZdvJVKtxItoaO7tB2nnAZnX2bzva14yIiBiv2nhsMa18wQ+ZbXv20+mlXwFr1Un3mdoN2643s3bxItsLJL0YOE/SNbZvaZH191P9oD/c9m2S1ge+3yIN0H7DyoiIiGileRCx0HbDno22t2t0TNK9kta2fXdpFnBfg9dYUP7eKukCqh/zTYMI29cDH6/Zvg04slmaIe02rIyIiIgWetgmYg6wd1nfG/j5s64trTY0bIKkaVRTebfs5ifpNkm3Dl/ayVRqIiIiIrrBPR0P4giqibI+CNwBvAtA0ixgP9sfAl4OfEvSEqpKgiNKLUMrtbUjKwC7A6u3k6kEEREREd3So4aVpUfkG+vsnwt8qKz/Afi7Ub52raMlXQ58rlXaBBERERFdMk5HpmxK0uY1m1Ooaibaig8SRERERHTJRJzFE/hqzfoiyoiY7SRMEBEREdENPRz2updsv2G0aRNEREREdIGYWDURkj7R7Ljtr7V6jQQRERERXaIlEyiKgGYDSrZVkAQRERER3dDbLp5dZ/vfASSdABwwNGunpNV4ZjuJhhJEREREdMuEqoh42sa1037bflDSs6atqCcjVkZERHRJD0es7KUppfYBAEmrky6eERERY6j1BFzj1VeBiyT9uGzvDhzeTsIEEREREV0gxn2NQ122Tyyziw7Nsr1bm8NlJ4iIiIjoGk/MqogSNLQVONRKEBEREdENE6x3RjekYWVERESXaHHjpaPXlXaXdJ2kJWXmzkbn7SDpJkk3Szq4s6u2liAiIiKiS+TGS4euBXYDLmx4bWkq8A1gR2AmsKekmR1fuYk8zoiIiOgG927ESts3AEhqdtoWwM22by3nngLswijaOrQrNRERERHd4iYLTJM0t2bZt8tXXwe4s2Z7ftnXM6mJiIiI6ALZrWoiFtpu1p7hV8BadQ59xvbPO81fLySIiIiI6JJO2j7Y3q7Dyy8A1q3ZnlH29UweZ0RERHRJn4e9vgzYSNL6kpYD9gDm9PKCCSIiIiK6wcBiN146IOntkuYDfw/8QtIvy/4XSDoTwPYi4GPAL4EbgFNtX9fRhVvo2eMMSStQdUVZvlznNNuHSlofOAV4PnA58F7bT/YqHxEREWOlV3Nn2D4dOL3O/ruAnWq2zwTO7E0unq2XNRFPANva3gTYFNhB0lbAkcBRtjcEHgQ+2MM8REREjBktccNlEPUsiHDlsbK5bFlMNcHHaWX/CcCuvcpDRETEmGnWvXMwY4je9s4oo2ddDmxINYrWLcBD5bkNNOnDWvrP7guw0lor8dqtejZWRtTx24t7OshZBAAbHHhxv7Mw6UzWe37bGFxDgDps+zDR9LRhpe3Ftjel6mayBfCyEaSdbXuW7VkrrLpCr7IYERHRNbIbLoNoTMaJsP2QpPOpWpWuKmmZUhvR8z6sERERY8KGAW370EjPaiIkrSFp1bK+IvAmqi4n5wPvLKftDYzLUbgiIiJGarI1rOxlTcTawAmlXcQUqv6qZ0i6HjhF0heBK4Hv9DAPERERY8NjNqjUuNGzIML2PGCzOvtvpWofERERMVgGtMahkcydERER0SWD2oCykQQRERER3TA07PUkkiAiIiKiC8TgduVsJBNwRUREdMuSJY2XDkjaXdJ1kpZImtXkvNslXSPpKklzO7poG1ITERER0Q0Getc741pgN+BbbZz7BtsLe5aTGgkiIiIiukQd1jg0YvsGAEk9ef3RyuOMiIiIbrBbPc6YJmluzbJvL3IBnCPp8h69/jOkJiIiIqJbmldELLTdrD3Dr4C16hz6jO12R3fe2vYCSWsC50q60faFbaYdsQQRERERXdLJ4wzb23V6fdsLyt/7JJ1ONbhjz4KIPM6IiIjoBlONWNlo6TFJK0laZWgd2J6qQWbPJIiIiIjoipZtIkZN0tslzaeaDfsXkn5Z9r9A0pnltOnA7yRdDVwK/ML22R1duIU8zoiIiOiWHg02Zft04PQ6++8CdirrtwKb9CQDDSSIiIiI6AYbFi/udy7GVIKIiIiIbjCweHLNBZ4gIiIiolsm2dwZCSIiIiK6wh03oJxoEkRERER0g0kQEREREaOUICIiIiJGbmwGlRpPEkRERER0g8Hp4hkREREjlnEiIiIiYtTSxTMiIiJGzpPucUbPJuCStK6k8yVdL+k6SQeU/Z+XtEDSVWXZqVd5iIiIGDM9nMVT0pcl3ShpnqTTJa3a4LwdJN0k6WZJB3d00Tb0chbPRcC/2p4JbAV8VNLMcuwo25uW5czGLxERETExmKphZaOlQ+cCr7S9MfBH4JDhJ0iaCnwD2BGYCexZ873bEz0LImzfbfuKsv4ocAOwTq+uFxER0Vc2eEnjpaOX9jm2F5XNi4EZdU7bArjZ9q22nwROAXbp6MIt9LIm4mmS1gM2Ay4puz5WqmSOl7TaWOQhIiKi13pYE1HrA8BZdfavA9xZsz2fHv94l3vcklTSysBvgMNt/1TSdGAhVc3PF4C1bX+gTrp9gX3L5kuBm3qa0camUeV3spmM5Z6MZYaUe7KZrOV+qe1VenkBSWdT3d9GVgD+VrM92/bsmvS/Ataqk+4ztn9ezvkMMAvYzcO+wCW9E9jB9ofK9nuBLW1/bDTlaUdPe2dIWhb4CXCS7Z8C2L635vi3gTPqpS03dna9Y2NJ0lzbs/qdj7E2Gcs9GcsMKXe/8zHWJnO5e30N2zt0mH67Zscl7QO8FXjj8ACiWACsW7M9o+zrmV72zhDwHeAG21+r2b92zWlvB67tVR4iIiIGgaQdgE8BO9v+S4PTLgM2krS+pOWAPYA5vcxXL2siXgO8F7hG0lVl36epWotuSvU443bgn3qYh4iIiEHwdWB54NzqNzoX295P0guA42zvZHuRpI8BvwSmAsfbvq6XmepZEGH7d4DqHJpoXTr7/kilTyZjuSdjmSHlnmxS7gnI9oYN9t8F7FSzfSZj+D3b84aVERERMZjGpItnREREDJ5JH0SUsSruk3Rtzb5NJV1chuWeK2mLsl+SjinDic6TtHn/cj56Dcq8iaSLJF0j6X8lPbfm2CGlzDdJenN/ct25JkOxry7pXEl/Kn9XK/sH5f1uVO7dy/YSSbOGpZnQ73mTMjccOniilxmalvsLpcxXSTqnPEcf+M94zfF/lWRJ08r2QJR7XLA9qRfgdcDmwLU1+84BdizrOwEX1KyfRdXWYyvgkn7nv4tlvgx4fVn/APCFsj4TuJqqQc/6wC3A1H6XYZTlXhvYvKyvQjV07EzgS8DBZf/BwJED9n43KvfLqcZguQCYVXP+hH/Pm5R5e2CZsv/Imvd6wpe5RbmfW3POx4Fjy/pAf8bL9rpUDQ3vAKYNUrnHwzLpayJsXwg8MHw3MPRL/HnAXWV9F+BEVy4GVh3WZXVCaFDmlwAXlvVzgXeU9V2AU2w/Yfs24GaqoVUnHDcein0X4IRy2gnArmV9UN7vuuW2fYPteoO4Tfj3vEmZGw0dPOHLDE3L/UjNaStR/R8HA/4ZL4ePouoaWdsAcCDKPR5M+iCigQOBL0u6E/gKSyc6GfMhRcfQdSwdY313lg5YMpBl1jOHYp9u++5y6B5gelkfuLLr2UPQ1zNQ5W5S5tqhgweqzPDscks6vPyfthfwuXLaQJdb0i7AAttXDztt4MrdLwki6tsf+Bfb6wL/QjVo1qD7APARSZdTVQc+2ef89Iyqodh/Ahw47Bcats0zf7EMjGblHlSNyqxq6OBFwEn9ylsv1Su37c+U/9NOAno2DHI/1Zab6v39NEsDpuiBBBH17Q38tKz/mKXVmmM+pOhYsX2j7e1tvwr4IdUzYRiwMqvOUOzAvUNVmeXvfWX/wJS9QbkbGYhyNyqzlg4dvFcJGmFAygxtvdcnsfRx5SCXewOq9i1XS7qdqmxXSFqLASp3vyWIqO8u4PVlfVvgT2V9DvC+0rJ3K+DhmmrwCU3SmuXvFOCzwLHl0BxgD0nLS1of2Ai4tD+57IxUfyh2qjLuXdb3Bn5es3/Cv99Nyt3IhH/PG5VZjYcOnvBlhqbl3qjmtF2AG8v6wH7GbV9je03b69lej+qRxea272FAyj0u9LtlZ78Xql/ddwNPUX3IPghsDVxO1Vr7EuBV5VwB36D6lX4NNS3aJ9LSoMwHULVo/iNwBGUgsnL+Z0qZb6L0WpmIS3lfDcwDrirLTsDzgV9TBYu/AlYfsPe7UbnfXt7/J4B7gV8OynvepMw3Uz0LH9p37KCUuUW5f0I1T9E84H+pGlsO/Gd82Dm3s7R3xkCUezwsGbEyIiIiRiWPMyIiImJUEkRERETEqCSIiIiIiFFJEBERERGjkiAiIiIiRiVBRMQYkfRYD15zZ0kHl/VdJc0cxWtcoGGzeEZEtCNBRMQEZnuO7SPK5q5UMzZGRIyJBBERY6yMkvdlSddKukbSu8v+bUqtwGmSbpR0UhmJD0k7lX2XSzpG0hll/z6Svi7pH4CdqSaOu0rSBrU1DJKmlaF/kbSipFMk3SDpdGDFmrxtL+kiSVdI+nGZiyAioq5l+p2BiEloN2BTYBNgGnCZpKFp2DcDXkE19PrvgddImgt8C3id7dsk/XD4C9r+g6Q5wBm2TwMo8Uc9+wN/sf1ySRsDV5Tzp1ENeb6d7cclHQR8AjisC2WOiAGUICJi7G0N/ND2YqrJv34DvBp4BLjU9nwASVcB6wGPAbfavq2k/yGwbwfXfx1wDIDteZLmlf1bUT0O+X0JQJYDLurgOhEx4BJERIwvT9SsL6azf6OLWPrIcoU2zhdwru09O7hmREwiaRMRMfZ+C7xb0lRJa1DVDDSbMfIm4MWS1ivb725w3qPAKjXbtwOvKuvvrNl/IfCPAJJeCWxc9l9M9fhkw3JsJUkvaadAETE5JYiIGHunU802eDVwHvApV9MT12X7r8BHgLMlXU4VLDxc59RTgE9KulLSBsBXgP0lXUnV9mLIN4GVJd1A1d7h8nKd/wP2AX5YHnFcBLysk4JGxGDLLJ4RE4CklW0/VnprfAP4k+2j+p2viJjcUhMRMTF8uDS0vA54HlVvjYiIvkpNRERERIxKaiIiIiJiVBJERERExKgkiIiIiIhRSRARERERo5IgIiIiIkYlQURERESMyv8HesaUO01MG/sAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "import s2spy.rgdr\n", + "\n", + "clusters = s2spy.rgdr._rgdr.masked_spherical_dbscan(\n", + " field_resampled.sel(i_interval=1),\n", + " eps_km=600,\n", + " alpha=0.05,\n", + " min_area_km2=3000**2\n", + ")\n", + "\n", + "clusters.cluster_labels.plot(cmap='viridis', size=3, aspect=3)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The above plot showing clusters will probably be useful for users to understand what is going on,\n", + "and will be required to finetune the parameters `eps_km` and `min_area_km2`. \n", + "\n", + "(Note: negative labels denote clusters with a negative correlation, positive labels clusters\n", + "with a positive correlation, and `0` labelled data is not within a cluster.)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`_rgdr.reduce_to_clusters` will both perform the clustering and return means per cluster:\n", + "\n", + "(and thus not providing direct feedback on what geographical region the clusters represent)" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
<xarray.Dataset>\n",
+       "Dimensions:         (anchor_year: 39, cluster_labels: 3)\n",
+       "Coordinates:\n",
+       "    index           (anchor_year) int64 1 13 25 37 49 61 ... 409 421 433 445 457\n",
+       "    interval        (anchor_year) object (1980-07-02, 1980-08-01] ... (2018-0...\n",
+       "  * anchor_year     (anchor_year) int64 1980 1981 1982 1983 ... 2016 2017 2018\n",
+       "    i_interval      int64 1\n",
+       "    target          bool False\n",
+       "    tfreq           int64 5\n",
+       "    n_clusters      int64 6\n",
+       "    cluster         int64 3\n",
+       "  * cluster_labels  (cluster_labels) float64 -1.0 0.0 1.0\n",
+       "Data variables:\n",
+       "    sst             (cluster_labels, anchor_year) float64 290.8 291.0 ... 298.2\n",
+       "    corr            (cluster_labels) float64 -0.3901 -0.06555 0.3602\n",
+       "    p_val           (cluster_labels) float64 0.01541 0.4659 0.02706\n",
+       "    area            (cluster_labels) float64 3.961e+06 3.869e+06 4.301e+06
" + ], + "text/plain": [ + "\n", + "Dimensions: (anchor_year: 39, cluster_labels: 3)\n", + "Coordinates:\n", + " index (anchor_year) int64 1 13 25 37 49 61 ... 409 421 433 445 457\n", + " interval (anchor_year) object (1980-07-02, 1980-08-01] ... (2018-0...\n", + " * anchor_year (anchor_year) int64 1980 1981 1982 1983 ... 2016 2017 2018\n", + " i_interval int64 1\n", + " target bool False\n", + " tfreq int64 5\n", + " n_clusters int64 6\n", + " cluster int64 3\n", + " * cluster_labels (cluster_labels) float64 -1.0 0.0 1.0\n", + "Data variables:\n", + " sst (cluster_labels, anchor_year) float64 290.8 291.0 ... 298.2\n", + " corr (cluster_labels) float64 -0.3901 -0.06555 0.3602\n", + " p_val (cluster_labels) float64 0.01541 0.4659 0.02706\n", + " area (cluster_labels) float64 3.961e+06 3.869e+06 4.301e+06" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "clustered_data = s2spy.rgdr._rgdr.reduce_to_clusters(\n", + " field_resampled.sel(i_interval=1),\n", + " eps_km=600, \n", + " alpha=0.04\n", + ")\n", + "clustered_data" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEXCAYAAACzhgONAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABs20lEQVR4nO2dd3hc5ZW43zPqvTerueKGsbENpncILQECJKSShDTSd5Nskt1N+SWbTTbZlE0PKaSQQgi992bAgG3cu2XZVu+9a77fH9+90mg05Y40TeZ7n2eembltjq5m7rmni1IKg8FgMBj84Yq1AAaDwWCIb4yiMBgMBkNAjKIwGAwGQ0CMojAYDAZDQIyiMBgMBkNAjKIwGAwGQ0CMoggTIrJbRC4Iss2/i8hvoyPR7BCR50Tkw7GWw+CbWP1/RORWEWkWkT4RKYj25881ROQCEamLtRyzxSiKMKGUWqmUei7INv+tlHL04xaRb4jIHWERLk4QkQtF5FkR6RaR2hD3/YCIbIyQaDH9fBGpFZFB6+JrP+ZF4rNClEuJyGKP90nAD4HLlFKZSqn22EnnH+t/Ne51Pi+ItVxzGaMoTlBEJDHWMvigH/g98MVof3Ccng9P3mpdfO1HQ6wF8kEJkArs9rUyzs7xK17n87lYCzSXMYoiTFh3hZcE2WbCShCR+dYd280ickxE2kTkP6x1lwP/DrzTuhvabi3PEZHfiUijiNSLyH+JSIK17gMi8pKI/EhE2oFviUiXiJzs8flF1p1rsYjkichDItIqIp3W64oInR4AlFKvKaX+DNSEsp+ILAd+BZxpnY8ua/lVIvKGiPSIyHER+YbHPvb5vUVEjgHPiEiCiPzAOtdHRORT1jaJ1j4+z6+/z48kofx/RGSxiDxvWWptInKnx7plIvKkiHSIyH4ReYfDz3/Bernd+pu/BOy3lnWJyDPWdkpEPikiB4GD1rKrRWSb9f17WURO8TjuqSKyVUR6ReROEfm7iPxX6GcoMtiWo4j8r3Xej4jIFQ72yxeR20WkwdrvviiIGzWMoog95wBLgYuBr4nIcqXUY8B/A3dad0OrrW3/AIwBi4FTgcsAT1fWBvRFuAT4JnAP8C6P9e8AnldKtaD/97cD1UAVMAj8zInAIvJu6yLg71EV+mnwj1JqL/BxJu8Sc61V/cD7gVzgKuBWEbnWa/fzgeXAW4CPAFcAa4C1gPe2f8DH+Q3w+VMQkV8EOCc7QvyzQ/n/fAt4AsgDKoCfWvJkAE8CfwWKgZuAX4jIimAfrpQ6z3q52vqb/wdYaS3LVUpd5LH5tejv3goRORVtNX4MKAB+DTwgIikikgzcB/wZyAfuAq73J4OInBPke3ZOgD/hVEtpHhCRr0po1s4GtFIsBL4H/E5EJMg+fwbS0eeoGPhRCJ8X/yilzCMMD6AWuCTINt8A7rBezwcUUOGx/jXgJu9trfclwDCQ5rHsXcCz1usPAMe8Pu8S4LDH+5eA9/uRbQ3Q6fH+OfRFMhLn6hKgNsR9PgBsDLLNj4EfeZ3fhR7rnwE+5iWHAhIdnt+Anz/L704f0GU97gvl/wP8CbjN87tkLX8n8KLXsl8DX3colwIWe7y3z2mi1zYXebz/JfAtr+PsRyvs84AGQDzWvQz8V5jP50JgAVrZrgL2AF8J4Xt2yON9uvU3lgbYpwxwA3k+1l0A1EXiexPNRzz5FN+sNHm8HgAy/WxXDSQBjR43Ny7guMc2x732eRZIF5ENQDP6YnMvgIiko+96LkffiQJkiUiCUmp8Rn9JlLH+ru8CJwPJQAr6LtUTz3MyD//ny8n5jSTXKqWest+E+P/5N7RV8ZqIdAI/UEr9Hv03bfBylSWi737Difd5vFlEPu2xLBl97hVQr6wrqMXRMMuCUsrTtblTRL6Jjot9x+EhJn6TSqkB6/vg73cJUAl0KKU6Q5V1rmAURfzi3db3OPqOt1ApNeZkH6XUuIj8A31n3Aw8pJTqtVZ/Hu3y2qCUahKRNcAbQDATGxF5D/rO1B8rlFLHgh0nRHy1Of4r2h1zhVJqSER+jHYX+NuvEe2asan0eB3s/AZtsywivwLe62f1UaXUSj/rfOH4/6OUakK71bDcMU9ZMYbjaFfjpSF87kzwPDfHgW8rpb7tvZGInA+Ui4h4KIsq4LCvg4rIucCjAT73CqXUiw7lC/q9ngXHgXwRyVVKdUXwc2KGiVHEL83AfBFxASilGtF+6B+ISLaIuERkkfXjC8Rf0S6I91ivbbLQfu8uEckHvu5UMKXUX9TUjBLvh08lYcmcir5zFxFJtfzW9vrnxCMg7UUzUOG5vfU3dFhK4nTg3UFE/wfwWREpF5Fc4Esef1Ow8+vr873Py8cDnJNQlIT9tzn6/4jIjTIZ6O5EXxjdwEPASSLyPhFJsh6niQ7O24Hb2gAyNKPdOKHwG+DjIrJBNBmikw6ygFfQMaDPWLK8HTjd34GUUi8G+Z75VBIicoWIlFivlwFfBe73WB/oexYy1nfnUXT8J8/6284Ltt9cwiiK+MV2obSLyFbr9fvRZvwe9AXhn2j/qF+UUq+ig77zmHp39mMgDWgDNgGPhUvwAJyHvvg9wmSA9gmP9ZXoOIovnkGnZTaJSJu17BPAN0WkF/gaWhEE4jfW5+1A350/gr5w2a6cQOfX1+dHkh/j/P9zGvCqiPQBDwCfVUrVWNbjZeggdgPapfI/aBcdBD7foONkf7QCx46ypZRSm9HWzc/Q5/AQ2u+PUmoEeLv1vgN9A3OPk+OGyMXADhHpR/+P70Enh9gE+7tnwvuAUWAf0AJ8ztdGIvKoiPy7x/s+y3JCRM61/odxh0x1FxoMscG6I/6HUuqsKH7mFcCvlFLV0frMeEJEnkArlb0xlOEP6GDvf0bp86L+PTsRMDEKQ1yglKoDIvrjFZE04EK0VVGCdufcG8nPjGeUUpfFWoZoE43v2YmIcT2FGcu07PPx+PfgexsijAD/D+0SeQPYi3ZZGQxB8fO7nnAdncgY15PBYDAYAmIsCoPBYDAE5ISLURQWFqr58+fHWgyDwWCYU2zZsqVNKVXka90Jpyjmz5/P5s2bYy2GwWAwzClExG+VvHE9GQwGgyEgRlEYDAaDISBGURgMBoMhIEZRGAwGgyEgRlEYDAaDISBGURgMBoMhIEZRGAwGgyEgRlEYDAZDNGg7qB9zkJgpChGpFJFnRWSPiOwWkc/62EZE5CcickhEdojI2ljIajAYDLPm/k/Bg9Muc3OCWFZmjwGfV0pttaZfbRGRJ5VSezy2uQJYYj02oAe3b4i+qAaDwTBL2g9BckaspZgRMbMolFKNSqmt1utedMvncq/NrgH+pDSbgFwRCTjRzWAwGOKO4V4YaIO+FpiDHbvjIkYhIvOBU4FXvVaVoweX29QxXZkgIh8Vkc0isrm1tTVichoMBsOM6LTaKI0NaqUxx4i5ohCRTOBu4HNKqZ6ZHEMpdZtSar1San1Rkc/mhwaDwRA7Oo9Mvu5riZ0cMySmikJEktBK4i9KKV9D1uvRg9BtKqxlBoPBMHforJ183dcUMzFmSiyzngT4HbBXKfVDP5s9ALzfyn46A+hWSjVGTUiDwWAIB1MURXPMxJgpscx6Oht4H7BTRLZZy/4dqAJQSv0KeAS4EjgEDAAfjL6YBoPBMEs6jkBuFXQdm5Oup5gpCqXURvSw+0DbKOCT0ZHIYDAYIkRnLZStgZ7GOWlRxDyYbTAYDCc07nFtSeQvhMxi6DWKwmAwGAye9NSDexTy5kNmibEoDAaDweCFHcjOX2ApirkXozCKwmAwGCKJrSjy5mvXk7EoDAaDwTCFzlqQBMiu0BbFQJuOW8whjKIwGAyGSNJxBHIrISERskpAuaF/brUaMorCYDAYIklnLeQt0K8zS/TzHHM/GUVhMBgMkaSzVscnwENRzK2AtlEUBoPBECmGumGww0NRFOtnY1EYDAaDAZiaGguTFkXv3GoMaBSFwWAwRArP1FiApDRIyTGuJ4PBYDBYeCsKmJO1FEZRGAwGQ6ToOAJpeZCaM7lsDlZnG0VhMBgMkcIzNdYmq2TODS8yisJgMBgihWdqrI2xKAwGg8EAwPgYdB/3oSiKYaQPhvtiItZMMIrCYDAYIkFPHbjHJlNjbewU2f65Y1UYRWEwGAyRwFfGE0wW3c2hAUZGURgMBkMk8KsoSvXzHEqRNYrCYDAYIkHHEXAlQXb51OVzsN+TURQGg8EQCTprIbcKXAlTl6fn6/kUxqJwhoj8XkRaRGSXn/UXiEi3iGyzHl+LtowGg8EwI3ylxoJWHBlFc6qWItYWxR+Ay4Ns86JSao31+GYUZDIYDIbZ03nEt6IAq+jOuJ4coZR6AeiIpQwGg8EQdgY7dYtx79RYm8wS43oKM2eKyHYReVREVvraQEQ+KiKbRWRza+vcGjFoMBhOQPxlPNlkFhuLIoxsBaqVUquBnwL3+dpIKXWbUmq9Ump9UVFRNOUzGAyG6QRVFJbryT0eLYlmRVwrCqVUj1Kqz3r9CJAkIoUxFstgMBgC03FEP/tVFKWgxmFgbnje41pRiEipiIj1+nS0vO2xlcpgMBiC0FkL6YWQkuV7/RwbiZoYyw8Xkb8BFwCFIlIHfB1IAlBK/Qq4AbhVRMaAQeAmpZSKkbgGg8HgDH+psTYTRXfNwMlREGh2xFRRKKXeFWT9z4CfRUkcg8FgCA+dR6DidP/r55hFEdeuJ4PBYJhzjI9Cd53/1FjwsijiH6MoDAaDIZx0HwflDux6SsmE5Mw5kyJrFIXBYDCEk2CpsTaZxcaiMBgMhjclwVJjbTJL5sxMCqMoDAaDIZx01kJCMmTNC7zdHGrjYRSFwWAwhJPOWsitBleQy2vm3GkMaBSFwWAwhJNAXWM9ySyG4W4YHYy4SLPFKAqDwWAIF0pB59HAqbE2cyhF1igKg8FgCBeDnTDc48yiyLJnZ8e/+8koCoPBYAgXnQ4znmBOVWcbRWEwGAzhwmlqLBjXk8FgMLwpcVpsB7q7LDInaimMojAYDIZw0VkLGcWQnBF824REyCgyFoXBYDC8qQjWXtybOVJLYRSFwWAwhIuQFcXc6PdkFIXBYDCEg7GR4O3FvZkjbTyMojAYDIZw0H0cUDOwKFrA7Y6UVGHBKAqDwWAIB6GkxtpklYJ7FIa6IiFR2DCKwmAwGMLBRLFdKK6nuVF0ZxSFwWAwhIPOWkhMnSykc4K9bW9TREQKF0ZRGAwGQzhw2l7ck4nq7PhOkY2pohCR34tIi4js8rNeROQnInJIRHaIyNpoy2gwvGloPQA9jbGWYm4yMgC1L0LZ6tD2myNtPGJtUfwBuDzA+iuAJdbjo8AvoyCTwfDm5C83wBP/EWsp5ia77oahblj3gdD2S8mCxLS4VxSJsfxwpdQLIjI/wCbXAH9SSilgk4jkikiZUsrc9hgM4WSwE7qOOms9YZjO67+FouVQfVZo+4nMiaK7WFsUwSgHjnu8r7OWTUFEPioim0Vkc2tra9SEMxhOGFr26uf2w3Gf0x931G+Bxm1w2i36wh8qc6DoLt4VhSOUUrcppdYrpdYXFRXFWhyDYe7RvFs/jw9bhWMGx7z+e0jKgFPeObP9s+K/31O8K4p6oNLjfYW1zGAwhJOWPZOv2w/FTo65xmAn7PonnHIjpGbP7BjGopg1DwDvt7KfzgC6TXzCYIgAzXsgf6F+bRSFc7b9FcaGYP0tMz9GZolWOGPD4ZMrzMQ0mC0ifwMuAApFpA74OpAEoJT6FfAIcCVwCBgAPhgbSQ2GExildIxi1Q3Q12oUhVOUgs2/h4rToeyUmR9nojq7BXIrA28bI2Kd9fSuIOsV8MkoiWMwvDnproPhbihZAYWLoe1grCWaGxx5XivV6349u+NklurnOFYU8e56MhgMkcaOTxSvhILFOvPJEJzXfwdp+bDi2tkdZw70ezKKwmB4s2NnPBUvh4IlOutpdDC2MsU7PQ2w72E49T2QlDq7Y82B6uyYup4MBkMc0LIHsisgLRcKFgEKOmqgZGWsJQuJL9+9g5313Zx3UhHnn1TEuuo8khIidC+89U+gxmH9h2Z/rAwrpd8oCoPBELc079HxCYDCJfq57eCcUhQtvUP8Y/NxynLS+M0LNfzyucNkpiRy1qICzl+qFUdFXnp4Pmx8DLb8ERZdPJkpNhsSkyG9wCgKg8EQp4yPQtsBWHKpfp+/SD/Pscynh3c04lbwxw+dRkl2Ki8fbuf5A608v7+VJ/boC/Ciogz+5dKTuPqUebP7sAOPQm8DXPW/YZDcIjO+i+6MojAY3sy0HdQT1mzrISUTsubNOUVx/7YGVpRls7g4C4C3rCzlLStLUUpxuLWf5w+08rfXjvG1+3dz6YoSUhITZv5hr/8WssthyVvCJD1x3+/JBLMNE+yo66J3aDTWYhiiyUTG04rJZQWL5pSiONY+wLbjXbxtzXRLQURYXJzJLecs4GtXr6Cjf4THds1iSFD7Yah5DtZ9EBLCeJ+dWQK9RlEY4pyh0XFu+OUr3P5SbaxFiU+Ugk2/hLotsZYkvDTvBlciFJ40uazAqqVQKnZyhcAD23VXn7euDuxSOmdxIVX56fzl1WMz/7DNv9fna+37Z34MX9htPOL0nBtFYQCgoWuQkXE3xzoGYi1KfPLGHfDYl+HVE2wkSssenRKbmDy5rHAJDHXBQEfMxHKKUor7tzVw2vw8ynPTAm7rcgnvOr2K1450cKilN/QPGx3U34NlV+tGfuEks0Q3ZBzqDu9xw4RRFAYA6rt03nxzz1CMJYlDWvbCI1/UrztqYitLuPHMeLIpWKyf2+O/QntfUy8HW/p425pp0wd8cuP6CpIShL++OoMOubvu0Qr0tFn0dfJHnI9ENYrCAGiLAqCp2yiKKYwMwF0f1EHeZVefWFXLQz3QfWxqfAI8FEX8xyke2N5Agku48uRSR9sXZqbwlpWl/HPLcYZGx0P7sJ136XMz/9wZSBqEiersWcRPIohRFAYA6juNovDJY1+G1r26n0/VmXPGJeMIe1iRd71EbjW4kuK+55NSige2NXDO4kIKMlMc7/fuDVX0DI3x8I4QG1F310HJyTMbThQMY1EY5gJ1lkXROzxG//BYjKWJE3b+E7b+Ec75F1h8sVW1zIljVbTYrTu8LIqERMhfEPcWxdZjndR3DXKNj2ynQJy5sICFhRn89bUQg9oD7ZBRGNo+TsmK7zYeRlEYgEmLAqDJxCm0Mnjwc1C5AS78D73MLkbriE9F0dA1SE8o6c3NeyA5C3Krpq8rWBz3iuL+bQ2kJLq4bKUzt5ONiPDuDVVsOdrJvqYeZzu5x/XMiPSCGUjqgNRcSEg2isIQ3zR0D1Kcpc335je7+2lsGP75QXAlwPW/g4QkvTyvGsQVdwFtt1vxmxdqOO97z/I/j+5zvmPLHt0I0JcrpWCx/jvdIfrxI0l3vZ6XAYyNu3l4RyOXLC8hMyX0eobr11aQnOjir05TZQc6AAXp/i2KOzYd5e+hWik2InFdS2EUhYFxt6Kxa4j18/MAaHyzK4onvw6N2+HaX0ydD5CYAjkVceV66uwf4cN/2sy3H9mLAo609TvbUSldQ+Gd8WRTsBjGR6BrFjUH4ebP18HD/wrAS4fbae8fCVo74Y+8jGSuWlXGvVvrGRhx4GodaNPPGb4tirFxN99/fD8/fWYWVlhmMfTG5wBPoygMtPQOMeZWrK3SiuJN7Xra94iuldjwcVh21fT1+YvixvX0em0HV/7kRTYebOOb16zkshUlztObext1YL7YT+M/uzlgvCjF9sPQtn/CmntgWwNZqYlcsLRoxod894YqeofHeHB7Q/CN+y1F4cei2Hy0k+7BUeq7BicyCEOmdBXUb4WxkZntH0GMojBMfLEXFWeSnZr45q2l6DoO990KZavh0m/63qZgEbTXxLSC1u1W/PzZQ9x02yZSEl3c84mzeP+Z8ynJTqWlx+Hc5WardUcgiwLip5bi0FP6uaeeodFxHt/dxOUrS0lNmnnPpvXVeSwpznTmfhpo189+gtlP7510GW0+2jkzgU66HEZ64ehLM9s/ghhFYaDOCmRX5KZRmpMaO9fTgcd1CmKsuO9W7ZO/4XbtZvJF/iI9NtS+cESZ1t5hbr79Nb7/+H6uXFXGg58+h5PLcwAoyU51nrXmL+PJJqMIUnLiJ6B98An9PNjJ87uO0jc8xjUOi+z8ISK8Z0MV2+u62VUfpCLadj35CWY/tbeFcxYXkp6cwObaGaZPLzgfElP17yDOMIrCMFGVXZ6XRkl2amwsioEO+NtNsPFH0f9sgI4jUPsinP/FyTRYX9jzB2IQ0N5U086VP3mR14508J23r+InN60hKzVpYn1JtpWM4OT/17wHssogPd/3ehF9HuKhlmJ0EGo3QoYuSntp6w4KM1M4c9HsM5CuW1tBapIreKpsv3Vj4ENRHG7t40hbP29ZWcLaqjw2187QokhOhwXn6TbmcdbzKaaKQkQuF5H9InJIRL7sY/0HRKRVRLZZjw/HQs4TnfrOQfLSk0hPTqQsJzU2RXdHngflnhzLGW0OPKafl7818HYxqqVQSvGpv75BRnIC933ybN51ehXila1Ukq1HcjY7cT+17NYZT4EoXBIfMYrajTA2BGveDcDR2gNcfUoZCa7ZF77lpCVx9SnzuP+NevoCWWIDbZCaM5kB58FT1ryLi5eXsK46j31NPaGlKXty0lugszY+FLQHMVMUIpIA/By4AlgBvEtEfNnBdyql1liP30ZVyDcJDV2DzLMaqpVmp9LaN8zouDu6Qhx+Rj8374nN3dT+R6FwafCJZbl2imx0L6BH2vpp6xvm4+cvYnlZts9tJhVFEEU/PgatB/y7nWwKFkNPHYw4zKSKFAef1C6ZU94JQOF4u8+W4jPlPRuq6B8Z5/5t9f43Gmj3G8h+em8LK8qymZebxmnz83EreONY18yEsWdc2DcucUIsLYrTgUNKqRql1Ajwd+CaGMoTU9xuxQ+f2M+hlr6of3Z91+BE582SnFSU0r7wcHO0vd93KqJScPg53b55uDv6cYqhbh1AXHp58G0Tk3WBWpTvtLdaF5611Xl+t3Hseuqo0Z1Kg406tQPasa4bOfSk7q9kKfGlGT2cWpkbtsOvqcxleVk2f331GMrfTUp/m0+3U2f/CJuPdnDJCl1ZvaYqlwSXsGWmcYrcSt0mJM7iFLFUFOWAZwvHOmuZN9eLyA4R+aeIVPpYf0Lw6pEOfvLMIR5wkqoXRpRS1HcOUp6nFUVZjr4rDXeK7Oi4m6t+spHfvXhk+sr2w7o53crr9Ht7mE60OPQ0uMfgpCucbR+DFNktRzvJSk1kcVGm320yUxJJT04I7noKFsi2iYfmgO2HtaJachmtQ0K7yuL0/MFpbrfZYAe1dzf0sL3OT1DbT/uOZ/e34FZwyXIdP8lMSWR5WRavzzROAdr9dOwVXQkeJ8R7MPtBYL5S6hTgSeCPvjYSkY+KyGYR2dza2hpVAcPFvW/ou+hoV0V3D47SPzI+aVFY7otwxymOdwzQNzxGQ7ePHPOaZ/XzmZ/Uz9GOUxx4DNLyoOI0Z9sXLNLB7yi6yLYe7WRtVR6uAH55EaE0O5Xm3iD/u+Y92n1WtDTwdnY8pi2GisJOi11yCY/sbKRRFbA41WHbjRC4Zs080pMTuGPTUd8b+LEontrbTEl2CifPy5lYtr46nzeOd87cfXvS5aDG9Q1MnOBIUYjIAifLQqQe8LQQKqxlEyil2pVS9u3Rb4F1vg6klLpNKbVeKbW+qGjmBTixYnBknEd26vbC0S52m8h48ohRQPgVxeFW7efu6PdRTHT4Ge37L1sD2RXRtSjc4zr1csllzkdb5i+E4Z7JIqwI0z04yoGWXtYFcDvZFGenBL/ZaNmjraKkwIN+SM7Qs6FjaVEcfELLmr+QR3Y20pdcTOZQ+NtcZKUmccO6Cu7fVj/xm5hAKZ8WxfDYOC8caOOiZSVTFPhp8/MZGnWzp2GGCq18nY6HxJH7yalFcbePZf+c5We/DiwRkQUikgzcBDzguYGIlHm8fRuwd5afGZc8ubeZvuExCjKSo56aajcDtF1P+RnJJCe4wi5HTauOvXQOeGWDjI/CkRdh0YU6JbNkxWQxWDQ4/po28U9yEJ+wiXJzwG3Hu1AKR4qixJFFEaB1hzcFi2NXdGenxS65FKUUexp6cOWWQ0+AoPMs+Nj5i1AKbnve6/863APu0WnB7FdrOugbHuPSFcVTltutcF6faZzClaBvXA4+oRMP4oCAikJElonI9UCOiLzd4/EBIHU2H6yUGgM+BTyOVgD/UErtFpFvisjbrM0+IyK7RWQ78BngA7P5zHjl3q11zMtJ5bKVpdFXFNbdk531JCKU5KSEveiuxrIoOr0tirrNuhp10UX6ffEKaDugFUg0OPCoDqIvvtj5PlFOkd16tBOXwGoHAVxdBzPsPyg70q/TL/217vDG7iIbi0w0Oy128aU09wzTOzxGUl6lVuwj4R/ZW56bxvVrK/jb68dp8VS2/b6L7Z7a20xqkouzFk1VICXZqVTmp828ngJ0nGKoC+pem/kxwkgwi2IpcDWQC7zV47EW+MhsP1wp9YhS6iSl1CKl1LetZV9TSj1gvf6KUmqlUmq1UupCpVQIrTHnBq29w7xwsI1rTi1nXk4qnQOjoU/emgUNXYOkJrkoyJicmVyanRp2F1hNmx+LouZZ7S9fcJ5+X7JS371FK498/2NQfbbOkXdKbhVIQtQsiq3HOllamu2oS2pJdiojY266B/0o2pZ9gArNohjqhv42vv3wHv74cq1juWeNnRY7/+yJbMDMYqslek9kkj5uvWARY+NufuuZdOGjfYdSiqf3tnDukiKfbUROq85n89FO/wo7GIsu0jcwcZImG1BRKKXuV0p9ELhaKfVBj8dnlFIvR0nGE5oHtjcw7la8/dTyiUCy4349YaDeqqHwzCIpzUkLu2Vjxyi6Bkam/ngOPwPz1upgMkymbEYjTtFRoxvNLXWY7WSTkKRbjkchbXTcrXjjWBfrqnMdbW+nyPpV9E4znmys5oCtR3fz241H+PYjezneEf67eZ8celLfQCSlcbClF4Cicsuai5D7aX5hBm9bPY87Nh2djKf5sCj2NvZS3zU4ke3kzbr5ebT1DXO0fYbnKjVb38DESZzCaYziOhHJFpEkEXnaqpZ+b0Qle5Nw7xt1rCrPYUlJFiURSk0NRH3nZA2FTWm2dj3N+G7Ii66BETr6RyjOSmHMrei1K2AHu6B+i45P2BQs0XdS0ch82m/drYUSn7DJXxQV19OB5l76hsccxSfAQXV28x5ISoc8h7kolptt5/bNKAUCfDeUmRczxU6LXXwpAAdb+shNTyKnpFqvj5BFAfDJCxczMDLO7S9ZVoUPi+Lpvc2IwEXLSnwe47T5ujXKjOMUoL+Xrft0hl2McaooLlNK9aDdULXAYuCLkRLqzcKB5l521fdw3am6fKTUaWVtGPEstrOx3Rdd3m6iGWJbE3aQr6vfOu6RF3TbDjs+AbqgrfCk6CiKA49C0TI99jNU8hfqC1mEffdbj2k/97oqPz2ZvCjJCvIdatmt/2aXw59+bjXKlURL7W5Om5/HrRcs4uGdjbO7ADrh4JP6ecklABxq7mNJcSaSbVVk90SuKHNJSRZXnFzKH16q1S48Hw0Bn9rbzOqKXIqyfDePXFyUSU5aEltm2kkWdJwCJhoiHm7t42Bz78yPNwucKgq7wclVwF1KqSCtFg1OuGdrPQkumWhHEG1FMTQ6TlvfyDRFUZaj34fLsrEzntZV64tdx4Bl0tc8C8mZ0+sXildE3vU01A1HX56ZNQH6TnukD/pawiuXF1uOdlKYmUxlfpBUVotiy/XU4u9/17zHeXwCwJXAcHY1+YNHefvaCj523iLKclL55oN7cLsjqCQPPTmRFgtwqLWPxcVZOqU3vSCiFgVoq6J3eIw/v1KrXU+JaTpdGH1ut9d1c+kK39YEgMslrKvOm51CLVikLewDjzHuVnzoD6/zhbu2z/x4s8CponhARPah6xieFpEi4E06tCA8uN2K+7fVc/5JRRRm6h93dloiqUmuqDXla+iamhprU5pj+bnDJMfh1n6SEoRTKnTAuNNWFIef0a0ZvButlayA7uP6Yh4pDj2lq7FDjU/YRClF1i60c1qJnJqUQG56km8l39ei746dZjxZHFHzWOhq4spVZaQlJ/Bvly9lZ303974RmTiBZ1osQHvfMB39IywutqrSs+fpsagR5OTyHC5aVszvNh5htLd1qttpn745uNhPfMJm/fw8Drf2+64dcsrSy6F2I8/uOMzR9gFq2vrD5hIOBaeKYitwGbAe+BLwF+BfIyXUm4FNNe00dg9NuJ3ASk2NQMaRP7xTY20mqrPDaFFU5adTZCnEzv4R7bbprJ3qdrKxL2QtESyb2f8YpOU7r8b2psBqHhjBOEVb3zC17QOO4xM2JVmpvmMUtjsvBItiZMzNpu485kszOclaWV2zupzVFTl87/F9zsaIhoqdFrtkMj4BsGRCUZRH3KIA+NRFi+kcGKWxsW6K2+npvc1U5KWxtCQr4P52nGJ27qfLYXyEN56/D4DeobGwuYRDwami+KpS6hhwJnAJ8H/ADyMm1ZuAu7fWk5WSOM18DWlK2Sxp8KrKtinOSkUkfBZFTVs/C4syyUvXKbidA6Nw2Grb4RnItrEvZJGKU4yPTVZju2Y4IS2nSgfdI5j5tNW6wISsKHJSfbuebHdeCBbFc/tb2DNaQiJjuh8X2q3ytbeuoLlnmF89H4G//+CT2tVTfY5+ayuKEk9FEVmLAmBtVR7nLC6kt6OZ8TStKAZHxnnxYBuXLC8JauWtKs8hOcE180FGAJUbGE/Opqr1RU5foBXPsWhlnXngVFHYif1XAbcppR4GkgNsbwjA4Mg4j+1q5MpVZdNysCNRw+CP+s5BXAKlOVNrJ5MTXRRkpIQlVjI27uZoez+LijLJSk0kwSXaoqh5FnIqJxvPeZJTCSnZkYtTHH9VFzM56Rbrj4REyJsfUdfT1mNdJCXIxAQ7p5Rkpfj+DjXv0ZPrMp23ublnaz0dqVanHY+eT+uq87n6lDJue+HwzGdE++PQk7DgXEjS38tDzb1kpiROxPDIngeDHREpuvPmUxctJnu8m9pB/dkvHWpjeMzNJcv9xydsUpMSWFWRM7s4RUIS21PXc3HCG/zbZTpV+WgcK4p6Efk18E7gERFJCWFfgxdP7Gmif2Sc69ZOb5ZbmqMVRTT8kHVdg5Rmp5KUMP1fWTqT6uy6LdM6XtZ1DjI6rlhYlIHLJeSmJdHdPwA1L8DCC3TbDm9E9FCdSLXyOPAouJJgUQjV2L7IX6jnZztAKcUn/7qVP71S6/jwW492snJeTshzoUuyU2ntHWbcO9jcvNN5/QQ6rfnpfc2sONlqsebV8+nLVyzDreB7j4UxXdYrLRZ0IHtRcebkHXxOhX7ubQzf5/phw4J8ihJ6eb3FxciYm6f2NpOVkjhxdx+M9fPz2FnfPeMi2uaeIe7oWE6hdHOy6O/asfbozwdxerF/B7rVxluUUl1APiY9dsbcvbWe8tw0Tp8//ctWnJUS1tTUQNR3Dk6LT9iUZodYdDfUA79/C/zlHTA2Gbw7bGU8LSrSGSN5GclkdezScyd8xSdsilfoVM5IKMz9j8H8s3VR02zIX+Q4RXZHXTcP72jkR08ecHTRGBlzs72uK2S3E+iiO7fSQeAJWg9A4/bJCngHPLijkdFxxeWnr9SV6149nyry0vnIuQu4b1sDbxwLU0tsr7RYgINWauwEdopsFOaWyNgQqWqIY0Pp3L21jqf3tXDe0iKSE51dOtdX5zM6rtjhr315EP70Si3Pjp+CEhepNU9SnJUy8yK+WeDor1VKDSil7lFKHbTeNyqlnoisaCcmLT1DbDzYyrWnzvPZMro0ikV3Dd2D0zKeJuXw477wx7FNuvVG3Wvw5FcnFts9nhYW6h96XnoSC7pfBURbFP4oWamznsIdtGw/rC94TmdPBKJgEYz2Q29T0E3v3Hwcl+j4zD1bg/vX9zT2MDzmnqGi8FF09/pvISEZ1t7s+Dj3bK1jWWkWK+bl6DRNH11kb71gMUVZKXzroT3hsYK90mK7B0Zp6R32UhSWJR6FgLZdlZ2WV8J/P7yX1t5hLnXgdrKx/38zcT8Njozzl1ePcdryxUjlBjjwGFX56XHtejKEiQe2N+BWcN2pFT7XR6uWYtytaOwamhbI9pSjK5S+U7Uv6AvRug/Cq7+CXbrhcE1bH/kZyeRZvaTy0pNZPrAF5q2B9ADmu93KI9wBbbt3zmziEzb22NQgAe3BkXEe3NbAtWvKObk8m99trAlag7BlhoFs8JG1NtwH2/8GK651HJ+oae3jjWNdvH1tuXb5FCz2OZciMyWRL162lK3HunhwxyxdQRNpsZdNLDrUqgvMJgLZMGlRRCGgbVdln3PKUnqHx0hwCRcsdR7jyc9IZnFx5owyn+55o46ugVFuOWeBLr5r2sGq7H6OxatFYQgfd2+tZ3VFzmROuBeO5x7PkpbeIcbcyq/rKeQBRrUbdarpFd+Dyg1w/6ehdT+HW/tZWJgxsVlpyijLxvfBQh/ZTp4UL7cEDbOi2P8oFC3XgejZUuCsluKRnY30Do/xjtMq+fA5Cznc2s/zBwMP2Np6rJPy3LSJ/0MoTPsO7bhTt8o+/aOOj3HfG/W4BK5ZY929Fy6G3gatdLy4fl0FK+dl8z+P7ptdQ8uJtNhJt9OhidRYj1TUpDSd2hwVRaEtitUnLWLlvGzOWlRAbnpoeTzrq/PYXNsRUoGi2634/cYjnFyereMhlgV8jtpKU89QVBuHglEUUWVfUw97G3um1E54Y1fWNnVHNkXWew6FNyFVZw91a//3/HN0C44b/6B/zHe+j8aWVhYWTSqKU8Z3kogb5Sst1pO0PO1iCGdAe7BLj5gMhzUBesiSKyloLcWdm48zvyCdDQvyuXJVGaXZqb5Hwnqw9WhnwPnYgSjMTMYlVnW2UvDab6BsNVSsd7S/26245416zllSNKmoJuZnT/9bE1zCV69eQX3XIL/ynuUQCl5psaDjE6lJrumWb050aino1xaFK7OIv3/0DH71Xp+z0wKyfn4+PUNjE2m+Tnj+YCuHW/u55ZwF2qIrWgq51axvuZt5tFHXGV2rwiiKKHLv1noSXcJbV8/zu01KYgL5GckRj1HYxXYV/lxPoVRnH9ukezbNt37g2fPght+h2g/yxZFfTLEolvW/zoBKYaDYwQ8u3K087GrscMQnwFGK7JG2fl470sE7TqtEREhOdPH+s6rZeKiNfU2+J6A1dA3S2D3EuqrcGYmVmOCiMDNFxyiOvgSte7U14bC6+/XaDuo6B7neMyuvQKdm+pt2d8bCAt62eh4/feYQrx2ZQTqoUnDw8SlpsaBrKBYVZU6P52WXR7w6G5js85RRQFZqEhkOWr17c5rV42zzUefn5fcbj1CclcJVq6xrhQhc+k0yB47zeMqXGH3t9qjOCDGKIkq43Yr7ttVzwdIiCjJ9NxKz0UV30VEUQV1PTuSofVHHJzyrnBdeQOPaz/O2hFe4oPveicVVXa+xyb2czhEHF62SFdC6P3xDjA48pitsHd5ZO6JgUcAU2X9sPk6CS7hh7WRM6t2nV5GWlODXqpiMTzhLwfTFRIX/a7dp6+zk6x3ve8/WejKSE7hsRenkQjseE2B+9revO5mq/HQ+/bettPWFaBE379KV+suumrL4UItXxpNN9rzoxSgkAVJzZ3yIqvx0CjNTHA8y2t/Uy4sH27j5rPlTs6tWXkvPB55np3sByzd/Ff58LXQdm7FcoWAURZRo7h2iuWeY85cG7g8Dus13xC2KzkFy0/3fIWWlJpGZkujMojjyolYSXjOYXyl7P0+Nn8rS7d/VI0e7jpHdX8tG9yo6+x1c/IutIUbhmNm8/e+w90HdEsFPNfbouDv0RncBUmTHxt38c0sdFy4totgj1pCbnsyN6yu4f1vD1ElqFluOdpKWlMCyssAtIgJRkp3CWFc97H0ITn1f8PnYFkOj4zy8s3Gir9MEyel6YNPhZ/SccR9kpSbx83evpWtglM/9fdv0Oo5A7H0IEFh65cSi/uEx6rsGWeKrVUZ2uS66Gw1zsZ83/W365sKhNeYLEeG0+c4bBP5uYw2pSS7es6Fq2rrc8iV8zPU1Hqr8op4O+Ysz4fXfRdy6MIrCwu1W/OjJAzyys5Ga1r7QvuQOsC+483KCBydLc1IjHsxu8NFe3JuS7JTgimKwC5p26OZ+XhxuG+Dfxj+hf9T/uBl23gXAC+5Vk40BAwoQhlYeI/1w3yfg3o/pAUkXf93nZqPjbi7+wfN8/4n9oR2/YCGMDfos/np2fyutvcO8Y33ltHUfPHsBo243d2yafkf4xrFOTqnI8VkI6ZSS7FTO731IuwRPu8Xxfk/s0fPb377WR1beeV+E45vg2W/73X/FvGy+ec1KNh5q46fPhDClcN9DUHUGZE7eSE3W4PiyKKKUIjvQPqUh4ExZPz+fus7BoL+ntr5h7tvWwPVrK3wGzUWEivxM7nZdBp94RVvHD/8r/Olt2iKLEEZRWDT1DPHTZw7yib9s5aIfPM/JX3+ca37+El++ewd/eOkIr9a00z2LIjj7C+Iki6UkO5W2vhFGxtwz/rxg+JpD4Y1dJR4Q7/iEBzWt/eTmFyPv/LP+wT39TcYySjmkyp0pisKTtNk/0zhF8x647ULY9ld9kbv5QcjynQP/3P5WjnUMcMcrR0NrdGd3kfUR0L7z9eMUZqZw4bLpVuSCwgwuXlbCHZuOTslgGRwZZ3dDz4zSYj0py0zgOveTjC+5LKQMr3u21lGem8YGX5XHa9+v6zBe/IG2zvzwjvWVvH1tOf/39EFeDJLdBejBPM27YNnVUxYfbPbq8eRJjq0oIux+si2KWbK+2lmc4o5NRxkZc/Ohc/zPSKm2aylyq+B998Fb/w/q34BfnKUTFyJgXRhFYTEvN40937ycBz91Dt+74RTedXoV6UkJPLa7iW88uId33raJNd96gid2By+u8oV9wfXuq+SLiZGoPtwS4UApFbAq28ZRdXbti5CQ4rMLa01bHwuLMnXGzVU/AGB8/vmA1e8pGIkpehRnqJlPSsGWP8JvLtQtRd53L1z0nzr47Ie7Nh8nJdFF7/AY928L4S51opZiqqJo6Rni2f0tXL+u3K9l8OFzF9DRP8J9Hu26d9R1MeZWs1YUawdepEi66Vj+fsf7tPQO8cKBVq47tdxnMSgAV34fytfBvbfqam8fiAj/de3JLCnO5HN/3xbcKt33kH5e7qUoWvpIShCq89On72NbFJEOaA+0hcWiWDEvm7SkBB7d2cTexh6fNyNDo+PcsekoFy4t8m1FWVQXpFPXMai9HiKw7gPauqjaELEZ26GH8E9g7CZeqyomm7AppWjpHWZPYw8f/uNmdtR1c9nK0gBH8U1TzxBJCUK+gxzsUo/K2oo8Hz+SWdI9OEr/yDgVflJjJ+TISaHF6hmU4O/CYddPJE1VgONuRW3bABfaMZm174OEZBIrTke27qHDqXVWvALqNzvbFmC4Fx76F+3mWnA+vP03fq0Im7a+YZ7Z18IHz57Piwfb+PMrR7nJylIKSk6FDuR7WRR3b61n3K18up1sNizIZ+W8bH638QjvtD5vi9UK49Sq2SmKFcfv5Ii7hI78M3BaHvbANqsY1EcPsgkSU+Adf4Jfnw93vgc+/LTPVijpyYn84j1redvPXuLTf9vK3z5yBon+XGl7H4KSVdMsn0MtvSwszPS9X1aZfo60RTHQHhaLIinBxVmLCnh4ZyMP79RuyqKsFKrz06kuyKC6IJ3OgRHa+ka45ZyFAY9VVZDOyLib5p6hyZu93Ep47z16mNYs4in+iKlFISKXi8h+ETkkIl/2sT5FRO601r8qIvNjICMl2alcuLSYkqwZNMqzaO4eojgr1f+dmgd+i+666/wGEkOh3k97cW9Ks1MZdyv/GSwT8Ynpbqe6zgFGxt1TaihY/U4SChaQk5ZElxPXE+gK7a5jupdUMBp36AvYrru1BfG+e4MqCdDFZWNuxY3rK3nfmdXsaexh67EuZ/K5EvT8aY/qbKUUd20+zunz8wPeGYoIHz53AQdb+njhoE7D3Hq0i4WFGeRnzKI5c+MOctu2cMf4pTT3OneX3r+tgdWVuQFlBrRyvPEPWjned6tfV8fi4iy+8/ZVvF7b6T/209eiu/l6WROgM578FaaSnG4V3UUwRjE+pi3S9NlbFAA/f89aHvjU2fz0Xafyxbcs5YKTinC5hJcOtfHDJw9w+0u1LC/L5uzFgRVTdb7+TU3r+SQCKTNPgAhEzCwKEUkAfg5cCtQBr4vIA0opTz/DLUCnUmqxiNwE/A+6g21MKMtNo7F7ZlkWTT1DjtxO4NHvyVMpNe2C286HVe+A6345Ixls7GK7oK4nu+iue8h3bOXYKzo+sWB6INvu8eTropOfnux86pfdyqNlrzat/dHTAH+4So+rvPkh3fTPAUop/rmljtWVuZxUkkV5bhrfeWQfd2w66tz9U7BoiqJ4vbaTmrZ+PnGhjxbqXly1ah7feWQfv32xhvOWFLL1WCcX+YhphMTrv0ElpnHX0HmUO7yxGRodZ09jD7eev8jZZyw4Fy77Fjz+77DxR3Cu7zlm16wp57UjHfz6+RpOq87nEu/xofseBtS0+MTQ6DjHOgYmK8N9Eem5FINWPCEMrifQHotTKnI5pSJ32rqh0XGOdwxQmJkS1JKtLtBehmMd/Zy5aPbWjhNiaVGcDhxSStUopUaAvwPXeG1zDfBH6/U/gYvF6UzICFCakzrjYT7NPcOOFUVeehLJCa5Ji0IpeOSLulhs+191G4pZUO9nBKo3tgvMrxVVu1HHJ8qn1yXYGSsLfSiK3PQk591x7bbYwVp5PP4fMDYMH3jYsZIA2FXfw76mXm5cp7N8MlISuX5tOQ/vaJzafTUQ+Qu1onDr5IM7Xz9OZkoiV64K7qJMTnRx81na5fX47mY6+kdmF58Y7IQdd8Ep72AwIYtmh3Gu/U29jLsVJ5eH0FH3jE/o+oxnvgWHnva72VevXsHKedl8/q7tHPduaLfvIe1ysm8ILGpa+3ErP4Fsm0jXUlgNAcPhegpGalICS0qyJnqiBaIsJ5VEl0S1i2wsFUU5cNzjfZ21zOc2SqkxoBuY9l8TkY+KyGYR2dza6iDLYobMy0mloXsw5C6ZSimauocmB68EQUQozvYYHLTzLjj2Mlz5v1ByMjz4WRiY+TCUhq5BUpNcFAT5UpZY1dl+A9q1L0Ll6dPiE6Cn2uWmJ/l0oeRnhGBR5FZBclbggPbhZ2H3PfqutsDhHbHFXVt0ENuzWv69Z1QzMu7mH5sdtrEuWKR7FPU20DM0ysM7G3jr6nmkJzsz2N+zQRfg/ed9u4CZNQKcYNtfYWwQOf0jFGc5n5a4u0G79lbOC2FIkgi87ae6d9bdt/hNz0xNSuAX71mL2634fw96KPyhbqh5XlsTXvd/B1usZoDFAVwpORGuzrYaAobLoggXiQkuyvPSotpF9oTIelJK3aaUWq+UWl9U5LyzY6iU5qQxNBr6rIieoTEGR8cdKwrwmHQ31ANP/KeuAVh/C1z7C/0FfmxaSMcx9V064ymYcVaYkUKiS3ynyA526ZiAj/gE6O6jnq07PMlNT3Yeo7CHGPlLkR0bhke+oOMEZ3/O2TEthkbHuX9bA29ZWUpOWtLE8iUlWZyxMJ+/vHrUWT2NnfnUfpgHtzcwNOrmnaf5D2J7k5uezPXrymnrGyYrNZHFwWIE/nC7dXpk1ZlQuiqkepzdDd1kpyYGTXCYRnIG3HSHdkHe+V6/U+eqCzJ475nVPLu/dbLrwMEndUHl8rdN2/5QSx8JLmF+YYBkDnvSXaSK7gaiZ1GESlV+elS7yMZSUdQDnr+mCmuZz21EJBHIAdqjIp0P7GK5UAPaEzUUDl1P9rbNPcPwwvegr1lbEy6XTjU974u6I+jeh0KSw6a+M3gNBejZyCXZqTT7+nuPvQIov4ricGu/36BoXnoSHU4VBejCu2Y/Q4xe/qmu3L7yf31aNoF4ck8z3YOj3Lh+enHZ+86YT13nIM8faAl+ILuWoqOGf7x+nKUlWayuCOHOHPjQ2Tpv/tSqPEcJDz45/Ax0HoHTPgxYBZMOFcWuhh5WzMt2lunlTf5CePtvdRzt8a/43eyGdRWMuxX32unAex+EzBKfqdWHWvqozk8nJTHAdL9s6/8WqYD2hOspviwK0HGKo1GcdBdLRfE6sEREFohIMnAT8IDXNg8A9qSVG4BnVDRmhPqhdEJRhHYHM1FDEaJFkd59CLXpl7oFQ4VHE71zPw+lq+Chz010twwFJ8V2NiXZfjK9jrwIiak+4xO9Q6O09g77jE+AnnI3NOp23iq5eKWece1d/dx5FF74X31H6tGa2il3baljXk4qZy2afiG4bGUJRVkp/PmVo8EPlF0Oiam0H9/L9rruiQaAobCwKJN7Tn6Z/1g4i+6rr/1aX3itO3SnrqexcTf7GntCczt5c9JlcNanYMsf4OgrPjdZVJTJ2qpc/rmlDjUyoC2KpVfqGyAvDgbKeLKJ9FwK2/UUaG5KjKjOz6BnaGxWRcChEDNFYcUcPoUesboX+IdSareIfFNEbFv0d0CBiBwC/hWYub8lDNhZQqFaFPYdeSiKoiQrma/we0jKgEu+MXVlQhJc+0vt/nk0tIm0Q6PjtPWNOFYUZTl+iu5qX/RZPwEeU+2KfLue8qxaEkfV2eDRysPL/fTYl0FccPl3nB3Hg8buQV482Mr16yp81ogkJbh41+lVPHegNbiJ73JB3gJaa/eQnOAK2EbeL0qx9ujtLH39a37dNwFp2AYHn9DWRKI+v6U5qfQNj9E3HLjSvKatn+ExNyvnzXI07AVfgZxKXcfip5HjDesqOdjSx5HXH9HTAX2kxY6Mualt6w8cyIbIt/Hob9PNABOSgm4abaqszKejHdGxKmIao1BKPaKUOkkptUgp9W1r2deUUg9Yr4eUUjcqpRYrpU5XSjmbZB8hCjO1z36mFoU9a8IJa/pe4JyE3bSe/kXfwbTSVXD+l3TNwO77HB+3IUjXWG/sLqRTDLnBTmja6bO/E+iKbJick+2NrSgcB7R9ZT7tewT2PwIXfEnn9YfIPVvrUUq7Q/zxrtMrcYnwl9eCWxXtKRUkdh3hspUlM6uBGOrWF87+Ftj8u9D3f+67+qK24WMTi0qygyQjWOxu0POcZ2VRgI5XXPl93db8lZ/73OTq1WWkJLro3HI3pOTA/OkzvI+29zPmVoED2RD52dlhKraLBHaKbLQyn06IYHa0SLB89iHHKHqGyEtPIjUpgL/Vk5F+Vu/+Hnvc1eyvuMH/dud8DsrW6KZgfc6yvZymxtqU5qQwMDJOz5DHXenRIPGJln4SXEJVvj9Foe/QHCcFpOfrSlzbohgZgEe/BEXLdIpmiEwUxC3Ip7rAt4ygralLl5fwj9ePB3STPbarkXuOplIlzXzt6mUhywNMuk+SMmDjj31OkvNL/RY48Cic9WlInbzYl2Q5m5a4u76HlESXX8UeEkuv0FlMz31Xuwa9yE5N4sqVhSzseJHxxZdOWD+e2FPtgrqektN1C/VIWRRhat8RCary7VoKoyjiktKcVBq7Qnc92cVrjnjxh6QMNPLV0Q/Q1BvAbWC7oIZ7tbJwEL6ZmGzn0KKw5Z5ysandqOMTfuY61LT1UZWfPrWXvgf2HbdjiwK0VWF3kX3xB9B9DK764YzcApuPdlLbPjBROxGI951ZTefAKI/s9D0P+h+vH+cTf9nKWO58khml2D3DXAv7Ynfhv+sL1Ou/cb7vs9/RVcoe1gQw0do8WJxid0MPy0qz/LfYCJUr/ke7BB/5os/v5Icqm8mjlzcyfN9oHGzpQ8RP11hvsisi6Hpqj8tANugWKYWZKVELaBtFESJlTjqqetHUM0SpU7dT+2F4+SeMr3onW9TS4OmNJSu0b3jvA7qWIAgNXYO4xFlzQpiMq0wpNKx9QddPJPr+m2q85mR7Y7dPdpwiC7ogq22/rtB+6f/glJtCKqzz5K7Nx0lPTuDKVWVBtz1rUQELizL486bpd8e/fbGGf7t7B2cvLuSDb71YLwwyP9svtvtkxTWw+FL9Nw73Bt/v+Gtw6Ek4+7PT2jfY/+NA3yGlFLsbulkxW7eTJzkVWuEdfHyy4Z8HK3teYJgkft3gu6fRwZY+KvLSps7D8Ef2POiJlOupDTLi0/UEduaTsSjikrKcVBq6Qiu6a3bavkMp7VJJSCHhsm+Sk5akU2SDcdZndEfPhz8Pvc0BN63rGqQkO9XxrINpimKgQ6dB+olPjLsVR9r6/QayQVdmA3Q4GV5kU7ISxkfgzvdBUrpuHzEDBkbGeHhHI1etKnM01lJEeO+Gat441sWueu3LV0rxv4/v578e3stVq8r47c3rSS05Se8QZH62X3oa9F14Vilc+BUdB3r118H3e/bb+q739I9MW5WZkkhGckLAG5u6zkF6hsZmH8j2ZsPHdaO/R780VeEphWv/I9TlncHTh/t8djo42NwbPD5hkz0vMhaFUlaMIj4tCtDtxqdVukcIoyhCpCwnjeEx50V3w2M6y8jJHAoOPKbvDi/4MmSVThbdBSMhcTILausfA27qtIbCxg7AT8gRpH6ioWuQ4TF3QLdBUoKLrNRE51lPMBnQbj8IF391yoCbUHhkZxP9I+PcGKCrqzfXr6sgNcnFHZuO4nYrvnr/Ln727CFuOq2Sn7zrVJ3rn1UGiWlTej6FRE+9Tm1NSNJK/6TLdY3IULf/fWpfgprn4Jx/0YFkH+ixuv5vNiYrssOsKBIS4eof6Yv4sx5ZaY3boPs4Oadeh1vB3VunWgPjbkVNW7/v8ae+yCnXF/TRMLfkH+rWLXPiNJgNOvOpsWeI4bHZNwoNhlEUIVJmWQYNDjOf7B9p0NRY97hO9yxaNuFrLgll0l3RUsirhtZ9ATdr6B50HMgG3X4hPyN5UlHY8YnydT63D9TjyZP8jOTQFEXRUnAl6uD9+g8538+LuzYfZ35B+sTAeyfkpCVx7Zpy7ttWz6f+tpU7Nh3jY+cv5DtvXzWZWuty6fPf5aDuwhc99ZPpnqDdiUNdsOlX/vd57jtauQQ4H1NawfhgT0M3LoFlpWFWFACVp+lZCa/+Ehq362V7HwJxUbjuWk6fn8/dW+qmWOfHOwYYGXOzyKmimEiRDXMtRZy27/CkuiAdpeB4R4THwWIURciU5U52VHWC/SMNWpXdU6975Wz42ESAttTJKFJPCk+CNt/DZEDfrTV2DTlOjbUpyfZohmj3dwoQnwD/NRQ2uenJdIZSLJSYAu/8i56F4GfmdTCOtvfz6pEOblhXEXJB3HvPqGZo1M0jO5v40uXL+MoVy6cfI6dy5sPuu+sn0z0B5q3R2UOv/Fxbit4ceUH/L875V53944fS7NSAjQF3N/SwuDjTWTxgJlzydX1X/tC/6JuhfQ9B9dmQUcAN6yuoaetnqzWDA3R8AnBuUUSqliKOq7Jt7KzCY1GopTCKIkQmLQpnF3DHVdl2QzW7HQT2SNRhxsYdjkQtPAnaDk10MfWmpXeIMbcKyfUEVgC/e8gjPjE9992mpq2P7NTEoA0H89KTnE2582Tp5fqufYb8c0sdIvieBx2Ek8tz+NSFi/nhO1Zz6wV+Gg/mVkHXcd/rAqGUvtB514Nc8GUY7oZNv5i+/bP/DVnz9B17AEqydSsYfzG13Q2zrMgORloevOW/dQrvE1/VFq/VUvzKVWWkJSXwzy2T7ie7GWDQ1FibiFkUlqKI82A2RKeWwiiKELGL7pocup6anFZl24rCY8pXSXYqbgVtfQ4vqIVLYGwQun1frCZSY0Ns/KYvNkNw9GUCxSdA11AsKs4Meseenx6i6ykM3PtGPecsLgzZorL5wluWBlYyuZW6SV0oNRCgXUyj/VMtCtBFlcvfBq/8Ymq34JpndazovM8H7W9VnJ3KiJ+YWlvfME09Q+GPT3iz6kY9bXCTVYS37CoAqxV7GQ9ub2RwRPvZDzX3UZaTSlaqw7TnSLXxmGjfEb+KoiAjmYzkBKMo4pGJojuHtRTNPUOkJLomMn380lmrffAefuqJjCOncYrCpfq57aDP1U4n23lTmp1Ke/8IY0de0AHb8rV+t61p62NhYfC7Qd1BNjp9akBfFOs6Bzn/pMh1FybHCpD7UdR+sd0mnjEKmwu+osdbvvIz/d62JnIqdQ+wIEyM1fXhfrID2SsirShErJqXZB1jyp1MJLhhXQV9w2M8bs2iP9TqoMeTJ5EqupsDricRoTI/PSpFd0ZRzICyHOfV2U3WwKKgPvHOWv3jT5hM2fQ56S4QhVaKpp84xUwVhe1uc9dsDBif6Bseo7lnOGh8AiA/I4m+4TFGxhy61WbJ/ibt0ohI0NYm13KLhRqnsGcq+FIUJStg5XU6Vba/HQ49BXWvw3lf8Pt/mLK7nbXm4zs00bqjLIKuJ5vCxfCuv8Fbfzxl8YYF+VTmp3HXluO43Srw+FN/ZEdgLsVAu07DDhD/iQeqC4yiiFtKc1Id93tq9jdG1JvO2mnD5f3OzvZHRoGu0G3zPZ+4vnOQ3PQkR/UDU+TISeU813aS23bD4ov9bndkYvxpcEUxo6K7WbC3Ud89LyuLzExhYPJOOVRFYbtNcvw0Ezz/SzDSDy//n66byK2GNe9xdOiSANXZuxt6qMhLIyeYtRsuFl8C806dssjlEq5fW8HLh9t5vbaDgZFx5zUUNpGYdNffFtfWhE11QQbHOgZwO5mbMguMonDK6BA8+TVo3c+83DQau4ccFd3pquyZKYqCjGQSXeJcUYBOI/XjemoIob24J/OSB/le0m30ZC2C0z/md7vJZoDB7wgn2nhESVHsb+qlMDOFwkznjRlDJqNYj4adiaIQF2T6GZ1avAxW3aDrKhregPP/zXHrkqIs/40B9zT0RD4+4YDr11agFHzvcX2DE7RrrDeRmJ0d51XZNlX56YyMuR2PvJ0pRlE45dEv6rYKD36OsuwUhsfcQdM7lVJaUQRLjR3u1aaul6JwuYTiLGfDZ0bH3bz7N5t4sjWH3ro9/Pr5wzy9t5natv6JKW32ZLtQqX716xTQw1PLvhUweHq4pQ+XTLZADoQds+kMpTp7Fuxr6mV5JK0J0LUUORUzi1Fklk5xO07j/C/p57wFun2JQ1KTEshLT5p2IekbHuNIW39kM54cUpmfzpkLC9hyVKfJhjzhLzsCRXdx3DnWk2hlPoXmg3iz8sZfYOufoPQUOPYyqypfBzJp7B4M2FK6c2CUkTF3cNeT3WXTR+qn06K7zbWdvHy4nQty53HpeCe/eHQz3egfXHKCiwWFGRxp6/c5pCcgu+4mee89/J96Bz1qAW8PsOnhtn4qg00ls7DPWzQyn8bG3Rxo7uX9Z848tdYxM0mR7a6bnvHkTeESuO42yF8QWKH4QNfBTHU92a64eLAoAG5cX8ErNe0UZiaTF2qbdttl19swOZZ2tvS3TyaHxDHVdi1F+wBnLIycYjMWRTCadunOrPPPhVuegLz5rNz7YwR30Myn2aTG2pR6FrsF4Jl9zSQnuHjfWy8F4KVbKrj71rP43vWn8MGz51ORl8aCwgzOXxpC1k9Po+4fVb6OB7JuCmrZBGsG6EnIw4uAbzywmxcOOGun7klt+wDDY+7IBrJtcmdQdNfT4D8+4ckpN/rt2BuI4uxUWrwsit31YZpBESYuP7mUzJTE0APZ4DGXIozupzhuMe7JvNxUElwS8QFGxqIIxFA3/ON9eiDMDb+HpDS48D9Iu+cjXO3aRGPPKQF3ty2B0pwgfvEAiqIkO5UXD7YFFfXpvS1sWJhPWpnugZTZW8O6tWezrtp5q4opKAUPfFqb89f9mqJ72gIqrFdr2jnc0sc5i53d1Uy6npwpir7hMf7wci1N3UOcF2KK676mKASybXKr9PCh0UH9fQmGUtq/vuTSiIlUkpXCfusc2Oxu6KEgI3kiKyrWpCcn8n83rQmeRu6LcFdnjwzA6MCccD0lJrgoz02LuOvJWBT+UAru/6R2C914+2QTupNvQJWs5AtJd9HcEbgN9ERVdrBZFJ21euBM2vSLekl28HGWNa191LT1c8nyEn2hSkgJ2MrDEVv+oBsUXvpNKFxCWU6aX0Xxj83Hee/vXqUiP40PnL3A0eFTEhPISE5w3Majtk3fMW073uVoe0/2NfaS4JKZ3a2GSk6VfnY6dW2oS1+UfKXGhonSnFRae4cnYlWgFcWKedkhtzKJJBcvL2Fd9QzmU08U3YWp3fgcKLbzpLog8l1kjaLwxys/h70PwqX/D6rPmlzuciEXf51qaab6WOD5D03dQ4hAcZYDi8KHNQGT1kigOMUz+1oAuGhZse6DVLB4doqiowYe/w9YeIGewYzVhbR3aEoa3rhb8Z1H9vJv/9zBhgUF3Hvr2SFlVeWmJzu2KGosRdHUMxRa/yu0RbGoKMNR7GTWhJoiO1FDESRGMQuKJyr8dZxiZMzNwZbeuHE7zZrkDG31h8uimGjfEf+uJ9CZT0eNoogBR1/RqbDLroYzPzV9/ZLL2Je0gotbbtdmqh+ae4YoyEgJPvshgKKYqKUIcHF8am8zS0uyqLTGI1IUuDlgQNzjcO+tukr8mp/rTB50g8LRcUW7dWHvHx7j43ds4dcv1PDeM6q4/YOnhZyPH0oHWbtGA0K3KvY29kYnPgHaogPnisK+uM1g7rdTSrxSZA809zI6ruImkB0WcsI46a7ftijmhqKoLkina2CU7sHIZRAaReFNXyv884M6A+naX+j2A96I8Gjpx8l3d8Br/ofL6NTYINaE261bU/uzKIK08egeGOX12k4uXu4xn6HwJK18xhwMPfLm5Z/C8U1w5femXLw8R6I2dA1yw69e4em9zXzjrSv41jUnOx6E5EluehIdDl1PR9r6KMpKISlBQlIUPUOj1HcNRic+AXouhSvReYqs7S6JoEUxOelOfx/2WK07Ti4/QSwK0OfPqbsvGHPOopjMfIoUMVEUIpIvIk+KyEHr2WfEVUTGRWSb9Xgg4oK5x+HuD+npYu/405RB9d4MlZ3Oc+5TURt/pLf3QVO3g2K73kY9uS2YReFn+MzzB1sZd6vpikK5Qx+i07RLV/4ufyuc8s4pq+yLzRO7m7jm5y9xvGOA33/gND5w9oIZ+7nzM5IdV2bXtPWzrDSL5WXZbDvu+3z7wm7dsTxaFoUrQccbQrEoAhXbhYESr5uNXQ3dZKYkUp0f3+0pQiKck+4m+jzNnRgFENHMp1hZFF8GnlZKLQGett77YlAptcZ6vC3iUj3737rP/1U/1J07A1Cak8r/jL4DGeqGl37ic5umHgftO+yMp1zfOf4ZKYlkpST6jVE8vbeZ/Ixk1lR66NrCJfq51XcrD58oBffdqpXj1T+eZknZCu8nzxwiNcnFPZ84iwuWzmzKnE2ewxiFUoojrf0sKMxgTWUuO+u6pwRmA7EvGq07vAmllsJJsd0sKchIxiXQYn2Hdjf0sLwsC5crfgLZsya7QlsC4Si6G2jXVmGAG8V4oio/8kV3sVIU1wD2zM4/AtfGSI5JWg/Aiz/QHTlPDd5Hpywnjb2qmq5F18KmX0Jv05T1Q6PjdA2MBrco7IlofiwK0EV3vgK4Y+NuntvfyoVLiycnrQEUWIrCTysPn7QdgKYdegaCD5O7KCuFwsxk1lfncd8nzuakktlfePPSk+kZGgs6b6Otb4Te4bEJRdE/Ms6hFmetvPc29ZKTluSsjUq4yK0KIZhd56yGYhYkJrgozNST7sbdir2NEZ5BEQts111vGKyKgTZtTcRRRlggMlISKcxMPvFcT0CJUqrRet0ElPjZLlVENovIJhG51t/BROSj1nabW1tDL8gCdAD4vf+EK7/vaHO7o+qupZ8C9yi8MHU/x5PtOmu16yGn0u8m/mZnbznaSffg6FS3E+iOlzlVoQW0azfq54UX+lyd4BKe/cIF/ONjZ1IQpn5JeRk6+N0VJAh3xMp4WlCYwerKXADH7qd9jT0sK82KbhpoTqV2KY45cKv1NEQ0PmFTmqMHGNW29zMwMh751uLRJieMtRT97XMmkG1TFeF24xFTFCLylIjs8vG4xnM7pTvr+fMjVCul1gPvBn4sIj5HiymlblNKrVdKrS8qmsW8gcWXOCuSAspytQI4Ml4Ea2/WdQceMQHbAihzoiiyKyDRf9uCicFBXjy9r4WkBOHcJT6+1IVL/HaR9UntRh2IDdACISs1Kazuionq7CDupyNWs8GFhZksKMggOzWRbce7gx7f7Vbsb+plWWkU3U5gZT6p4Hn9drFdduQynmyKs/R3yJ5BcUJlPMFkHUo4qrPnSENAT+wuspEiYopCKXWJUupkH4/7gWYRKQOwnlv8HKPeeq4BngNO9bVdLCjM0Bk4Dd1DupunKwme/c7E+pBGoAYZ71mSnUJL7/C0VsJP723mjIUFvqeBFZ6kXU9+xqJOQSk4+pKeZRzFO+/JNh6BLYqatn6SE1yU56XhcgmrK3MdZT7Vdw3SPzLOsrIoXxQnaimCxCkmiu0ib1GUZKdYiqKbpAQJvZV3vBPOSXcDc9OiaOgeZHhsPCLHj5Xr6QHgZuv1zcD93huISJ6IpFivC4GzgT1RkzAILmvSXVP3EGSVwhkfh5136cwhQnQ9BYhPgHYbjLsVbf2TmU+1bf0cbu3XRXa+KDpJX4Sc+GzbD0Nfc8ARp5HAdj11BLMoWvupLkifiMOsqcxlf1MPAyP+q9XBYwZFTCwKgqfI2ne/EY5RgLZKOwdGeeNYFyeVZJGceIJlxk8U3YVBUfS3zZmMJ5vqgnSUgrpOZ3NyQiVW35bvApeKyEHgEus9IrJeRH5rbbMc2Cwi24Fnge8qpeJGUYB2KzVYU+M4+7M6NXK3rtZu6h4mPTmBrEBDgkYG9AU6iKKYLLqbVBRP7W0G0G07fBFk2t0Ual/Uz9FWFA6HFx1p0xlPNmsqc3Er2FkX2P20r6kXEcISeA+J7HIddwoW0A40AjXM2JbtlqOdJ57bySanEhq3awt5poyPaktvjtRQ2NgpspEKaMdEUSil2pVSFyulllguqg5r+Wal1Iet1y8rpVYppVZbz7+LhayBKMtJmwwyp+VB8Qqo3wpAU88gpdlBRqA6yHgC30V3z+xr4aSSzMlqbG9sRdHqQFEcfUkP3SlYHHzbMGIrikDDi8bdiqPtAywomqooALbXdQU8/r6mHqrz00Oe6DdrEpJ0vCeY62mi2C7yiqLYav437lYnXsaTzdr36zGxu++d+TEGOvTzHLMo7KK7o+2RqaU4wezP6GLPzp6YdFe+Fhq2glI0ORmBGqBrrCeTlbVaUfQMjfLakQ4uWuYvWQzIKNKmeDCLQikdyJ5/TtTTAdOSE0hNctEVIEbR0DXIyLh7SvvygswUKvPTgsYp9kWzdYc3TlJkJ4rtAvwfw4Tnd/GEtShOu0XPjHn83/UwsJkwx6qybQozk0lPTuBYx4nlejohKMtJZWTMPeljL1+nW5N31NDcMxx8sp1DRWEXTNmK4vn9rYy5FZd4p8V6ImIFtIMoio4anco5/+zA20WIvPTkgDGKmonU2KmdX1dX5LLtWJff/QZHxjnS3h/dQjtPciqhO4ii6K7XlkcEi+1sbKtUBJZHO7gfLVwJcPWPdE3Tc9+d2THmWOdYGxGxUmSNRRF32P2PGu1iuHlrAXDXbabZUVX2UUjODPqlTExwUZSVMpFy+8y+FvLSkzi1KsisCTvzKRBHX9LP888NvF2EyEsP3MbjSKtOjV3gNRBpTWUuDd1DE9XG3hxo7kUpYmtRdNfDeICAe099VDKeQPfVSk5wsaAgI/quuGhSsR7W3ayLYK3EkpCYaN8xtywK0LPqR8dnEZ8JgFEUs2CeVUsxoSiKlkFSOkNHNzPmVs5qKPLmO3L52EV3Y+Nunt3fMr0a2xeFS6CvSVs5/qjdqN1UdkwjyuRlJAW0KI609ZNlVZ56cmpVLuC/k6w9rCjic7L9kVsJalxba/7oqY9KfAL0HeeCwgzWznSQ1Vzi4q/r9hsPf95ZergntkUxx1xPAD9796n88UOnR+TYRlHMAtu11Nht+QUTEqFsNdRtAXAWowjidrKxi+7eON5F18AoF/vLdvKkyJr568+qUApqX9LzNmLUrkBbFP5jFDVt/SwoypiWFLByXg6JLv+dZPc29pKenEBlXowa3wVrN66UVZUdHUUB8OcPn87X37oiap8XM9Lz9cCt45tg+99C29e2KNJmMEApxkSy+4BRFLPALrpr9OzDVL6OlPZdJDIWOEahVEiKwm7B8NTeZhJdwrknObjjCZYi21mrM29i5HYCqzFgINeTV2qsTWpSAsvKsgJaFEtLY9j4LidILcVgp65ziUINhU1xVqrv4swTkTXvgcoN8ORXJzOZnDDQpjMYoxA3mksYRTEL7KK7xi6PTIPytSSMD7NU6gJXZfe1wNhgSBZF9+Aoj+5sYsPCfLKd/OBzqyEh2X8XWTs+UR2bQDZAXkYyXYOjPrvBDo2OU9816FNRgI5T7KjrnlaxrpRiX1MMM55gcpaHP4tiooYiOjGKNx0uF1z1A62Qn/mW8/0G2udcIDsaGEUxS+blpE21KKyA9hrX4Wl+9Sk4zHiysd1YxzoGuDhQWqwnCYmQv8i/66n2Jf2jKFrm7HgRIC89SXthfDQGPNo+gFLTA9k2ayrz6Bse43Dr1E6yzT3DdA2Mxi4+AZCUqtNe/SoKewRq5Ps8vWkpXQUbPg6bb4f6Lc726W+bk4HsSGMUxSwptWopJsibT39CDqen1JIYaOpbiIrC0zqZ1i02EIVL/Lueajfq+IQrdl+DQEV3ns0AfWEX3r3h5X7a22S37ohxGmhOpX/X04SiMBZFRLngK1phP/SvejBZMAba52QgO9IYRTFLynJ1v6eJojsRDiUtYbUcDryjrSgCtBf3xB6purg4k+oC33fYPik8CTqP6NYEUz7/qM7zr45u2w5v8jL8t/GwayjmF/oOSC8szCArNXFanGJfoy62WhrtHk/eBCq6664HSdB9wgyRIzUb3vJtaNwGm38ffPs52OcpGhhFMUvKslMZGXfT7pHiucO9mKrxYzASoPilsxay5mkXhQNKc9JIcIn/3k7+KDwJ3GPTx6JO1E/EWFGk240Bp7uejrT2U5SV4jcA63IJqyty2e6tKJp6KM9NIyctxoHb3Eo9mMhXimZPg1YSroToy/Vm4+TrYcH58PS3dGzQH0qZGIUfjKKYJWW5uujOcwLdppFqXLh1gzJ/dB117HYCyExJ5M6PnsGnLwqxH1ORn8yn2pcm+1PFkMlW475cT/1TWnf4Yk1lLvuaehkcmXQrxGQGhS9yq/Q89L7m6et66qKaGvumRkQHtkcH4KF/8e+CGurStS/G9TQNoyhmiV1UZ3eRHRgZ49Wh+XploABaCKmxNuvn54deVTsxFtVLURzdqLOdYhifgEnXk6/hRUfa+llYFFxRjLsVuxp0UeHImJtDLX2xa93hSaAU2ShNtjNYFC6By74F+x6CBz/j28rrt9t3GEXhjVEUs6TMauNhd3Zt6h6ijRwG0somOslOY3RIXyhCVBQzIiVT37l6dpHtrtOKKoZpsTYZyQkkJ7imDS/qHhilvX/Eb8aTzcRoVKvv0+HWPsbcKvaBbPBfdKeUjlHkmIynqHLGrXD+l+CNO3TjQO925BMNAY3ryRtTVTJLCjKS9aS7LktRWApjoGgN6f4siu7jgIqOooDpmU+1dnwi9opCRMhNT5pmURxp990M0JuirBTKc9PYZrUc3zeR8RQHFsXEpDsvRTHYqWtojEURfS74Cgz3waaf65uoi/5zct0c7vMUaYxFMUtcLqE0J5Umq42H3eFVytfqOIRtznoSYmrsrClcqmsp7Duooxt1L5ySk6Pz+UHIz5henW2nxgazKADWVE12kt3X2Kub3znYL+IkZ+jAqLeiiOLAIoMXIjoLau374YXvw8YfT66bo51jo4FRFGGgLDtNz85mskFg+gKrOVeDD/dT1BXFEhjp1e2XQddPVJ0VNxk3uelJ0/o9HWntxyV6FnAw1lTkUt81SGvvMHubellSkhm4hiWa+KqlmKihMIoiJojA1T+Gk2+Ap74Or/1GL5+jsyiiQZz8muY2di0FQHP3EFkpiaRVrwPEd5yisxYS0yAzhMK52TDR82m/vpvtqIl5Wqwn+RnJ0wruatr6qcxPdzTbeY3VSXb78S72NfbER3zCJrfSh0URvVnZBj+4EuC6X8HSK+GRL8C2v2nrPykDktJiLV3cYRRFGNCupyHcbkVTzxAlOamQkqW7t/qKU3TWQl519Dq2enaRjaP4hE2uj5kU/poB+uLkeTkkuIRn9rfQ0jsc29Yd3uRW65GonoFTu9guCpPtDAFISIIbbtc1Fvd/Ag48agLZfjCKIgzMy0ljZNxNx8AITT3Dk3MoytdpReGdXTGD1NhZkVkCKdk6oH10o35dekr0Pj8IeelJdA6MTlS3K6VCUhRpyQksK83igW3a9x9XFkVOpQ5c24FSsIrtyuLG9femJikVbvorVJymLW0Tn/CJURRhYGIuRdcQzZ6zsuedqv2enj7qENuLhwWRycyn2peg6sy4ukjlpScz7lb0DOlpcC29wwyMjActtvNkdWUufcN6/7ioobCxU2Q9x6L21JmMp3giJRPe/Q8oX6/nyRimERNFISI3ishuEXGLyPoA210uIvtF5JCIfDmaMobCPKuWoq5zgNa+4ckGfuXr9LOn+2mgA0b6oqsoQMcp6rdC+8G4cjuBR3W2lSJb0+osNdYTu0FgYWYKhZkp4RVwNkykyHrcLPQ0mPhEvJGWC7c8qYPchmnEyqLYBbwdeMHfBiKSAPwcuAJYAbxLROJyPJdtUeys72bcrXSMAnT6aULy1IB2tDOebApPgmFdYxBPgWzQwWyYbONRY6fGBqnK9uRUS1HEVXwCJps+2gFtu9jOZDzFHy5XzCY9xjsxKbhTSu2FoKP7TgcOKaVqrG3/DlwD7Im4gCFSkJFMcoKLN6xc/gmLIjFZ98SfoiiO6OdYKAqA5CwojS/zOtdqDGgriiOt/aQkuigLNkrWg0VFmZRmp7Iu3mZCp+VCSs6k+9EU2xnmIPFcmV0OeCag1wEbfG0oIh8FPgpQVVUVecm8cLmEkpwUdljVwVMm281bC9v+qhuRuRImLYrc6ugKaSuKqg1xN+ZxwqKwOsjagexQxpi6XMKT/3oeqUnxE3uZwLPduKmhMMxBIuZ6EpGnRGSXj8c14f4spdRtSqn1Sqn1RUVF4T68I8py0ui3OpiW5Hj4yMvXwWj/ZAuNzlqdhZQcvJAsrOQv0Jk2y98a3c91QK5XB1knzQB9kZWaRFK8FNp5kls5GaMwVdmGOUjEbi2VUpfM8hD1gOdUnwprWVxip8QmuoTCDE9FoUejUr8FipdrRRFtawJ0zvi/7o3+5zogOzWRBJfQOTDC6LibYx0DXLHqBBrok1MJR1604hN11jKjKAxzhzi8/ZrgdWCJiCwQkWTgJuCBGMvkF7uLbEl26lSXScESHRew4xSdoc2hCCsicRmsExHy0pPo6B+lrnOQMbcKKeMp7smt0i1Uhrq0RWGK7QxzjFilx14nInXAmcDDIvK4tXyeiDwCoJQaAz4FPA7sBf6hlNodC3mdYFsUJdleqZkuF8xboy2KsRGdQx8rRRHH5FnV2aE0A5wzeHaR7ak3xXaGOUessp7uBe71sbwBuNLj/SPAI1EUbcbYisJOlZ1C+Tp45efQcRiU2ygKH+Sl6w6ydg1FKMV2cc/EXIrjWlGYjCfDHCOeXU9zCk/X0zTK14F7FPY+pN8bRTENPZNilCNt/eSmJ01MvjshyPEYYNRdb+IThjmHURRhojwvDZdAZZ6PbCY7oL3rbv1sFMU07JkUofR4mjOk5+uupN3HrRGoRlEY5hbxlVA/h8nPSObOj53Jynk+GtJll+vgZeteXamdVRZ9AeOcXA/X01mLT7DGbCI6TtG43Sq2M4rCMLcwFkUYOW1+PunJPnSviC68A50a6zKn3Zv8jCRGx3Wb9hMqPmGTWzXZ88vEKAxzDHPFihZ2g8C8GNRQzAHsojsIrRngnCGnEsaGrNcVsZXFYAgRoyiiRfmp+tnEJ3ySP0VRnKAWhY2xKAxzDKMoosW8tZCUrjvKGqaRl5E08Xp+YZTbm0QDu5bCFNsZ5iAmmB0t0vPhM9vMBC0/2DMpynJSfcd55jp2iqwptjPMQU7AX2Qck2XuJP1hK4oT0u0Ek64nU0NhmIMY15MhLshOS8IlJ7CiyCiChBQTnzDMSYxFYYgLElzC929YzZqq3FiLEhlcLjj383qOusEwxzCKwhA3XL/uBE8bveBLsZbAYJgRxvVkMBgMhoAYRWEwGAyGgBhFYTAYDIaAGEVhMBgMhoAYRWEwGAyGgBhFYTAYDIaAGEVhMBgMhoAYRWEwGAyGgIhSKtYyhBURaQWOzuIQhUBbmMSJBEa+2WHkmx1GvtkRz/JVK6WKfK044RTFbBGRzUqp9bGWwx9Gvtlh5JsdRr7ZEe/y+cO4ngwGg8EQEKMoDAaDwRAQoyimc1usBQiCkW92GPlmh5FvdsS7fD4xMQqDwWAwBMRYFAaDwWAIiFEUBoPBYAjICa8oROT3ItIiIrs8lq0WkVdEZKeIPCgi2dbyJBH5o7V8r4h8xWOfy0Vkv4gcEpEvx6F8tdbybSKyOUbyJYvI7dby7SJygcc+66zlh0TkJyIicSbfc9b/d5v1KA6TfJUi8qyI7BGR3SLyWWt5vog8KSIHrec8a7lY5+eQiOwQkbUex7rZ2v6giNwch/KNe5y/B2Ik3zLrfz8sIl/wOlbYf8Nhli8iv+GwoJQ6oR/AecBaYJfHsteB863XHwK+Zb1+N/B363U6UAvMBxKAw8BCIBnYDqyIF/ms97VAYYzP3yeB263XxcAWwGW9fw04AxDgUeCKOJPvOWB9BM5fGbDWep0FHABWAN8Dvmwt/zLwP9brK63zI9b5etVang/UWM951uu8eJHPWtcXB+evGDgN+DbwBY/jROQ3HC75rHW1ROA3HI7HCW9RKKVeADq8Fp8EvGC9fhK43t4cyBCRRCANGAF6gNOBQ0qpGqXUCPB34Jo4ki9ihCjfCuAZa78WoAtYLyJlQLZSapPSv4g/AdfGi3zhkCOAfI1Kqa3W615gL1CO/v780drsj0yej2uAPynNJiDXOn9vAZ5USnUopTqtv+vyOJIvIoQqn1KqRSn1OjDqdaiI/IbDKF9cc8IrCj/sZvJLciNQab3+J9APNALHgP9VSnWg//HHPfavs5bFi3yglcgTIrJFRD4aQdkCybcdeJuIJIrIAmCdta4cfc5sYnX+/Mlnc7tl9n81XK4xT0RkPnAq8CpQopRqtFY1ASXWa3/ftYh/B2cpH0CqiGwWkU0icm04ZQtBPn/Ey/kLRDR/wyHxZlUUHwI+ISJb0ObiiLX8dGAcmAcsAD4vIgvniHznKKXWAlcAnxSR82Ig3+/RP8DNwI+Bly15o81M5HuPUmoVcK71eF84BRKRTOBu4HNKqSlWoGVlxTRPPUzyVSvdnuLdwI9FZFGcyRcxwiRfNH/DIfGmVBRKqX1KqcuUUuuAv6F9l6C/4I8ppUYt18RLaNdEPVPvPCusZfEiH0qpeuu5BbgXrVSiKp9Sakwp9S9KqTVKqWuAXLTPth59zmxicv4CyOd5/nqBvxLG8yciSeiLyF+UUvdYi5ttl4313GIt9/ddi9h3MEzyeZ7DGnTM59QYyOePeDl/fonmbzhU3pSKQqyMFhFxAf8J/MpadQy4yFqXgQ7W7UMHR5eIyAIRSQZuAsKS1REO+UQkQ0SyPJZfBuzyPm6k5RORdOvzEZFLgTGl1B7LBO8RkTMsl877gfvjRT7LFVVoLU8CriZM58/6e38H7FVK/dBj1QOAnbl0M5Pn4wHg/VZ20RlAt3X+HgcuE5E8K4PmMmtZXMhnyZViHbMQOBvYEwP5/BGR33C45Iv2bzhkwhkZj8cH+o6yER08qgNuAT6LvpM8AHyXyQr1TOAutI97D/BFj+NcaW1/GPiPeJIPncmx3XrsjqF884H96IDeU2hXhH2c9egv/mHgZ/Y+8SAfkIHOgNphnb//AxLCJN85aLfDDmCb9bgSKACeBg5asuRb2wvwc+s87cQjEwvtUjtkPT4YT/IBZ1nvt1vPt8RIvlLre9CDTlaoQydSQAR+w+GSjwj+hsPxMC08DAaDwRCQN6XryWAwGAzOMYrCYDAYDAExisJgMBgMATGKwmAwGAwBMYrCYDAYDAExisJgMBgMATGKwmCYISIyXzzamxsMJypGURgMMcLqAnzCfZbhxMMoCsObFhG5z+rUudvu1ikifSLybdGDjTaJSIm1vERE7rWWbxeRs6zDJIjIb6xjPCEiadb2a6z9d1j72YNrnhORH4seTPNZHzJlicgRq5UIIpJtvxeRRSLymCXziyKyzNrmrSLyqoi8ISJPecj8DRH5s4i8BPw5wqfTcAJjFIXhzcyHlG4cuB74jIgUoNt5bFJKrUbPtPiIte1PgOet5WvRbRYAlgA/V0qtRLdksGdf/An4klLqFHRLi697fG6yUmq9UuoH3gIp3ZTwOeAqa9FNwD1KqVHgNuDTlsxfAH5hbbMROEMpdSp6zsK/eRxyBXCJUupdIZ0Zg8EDY44a3sx8RkSus15Xoi/6I8BD1rItwKXW64vQzQxRSo0D3ZaVcEQptc1j+/kikgPkKqWet5b/Ed2jy+bOIHL9Fn2xvw/4IPAR0W2szwLukslRGSnWcwVwp9WlNBk44nGsB5RSg0E+z2AIiFEUhjcloudlXwKcqZQaEJHngFRgVE02QBsn+G9k2OP1OHryYDD6A61USr1kBcovQDcn3CV67neXUmqNj11+CvxQKfWAtc83nH6WweAE43oyvFnJATotJbEM3bI9EE8DtwKISIJlNfhEKdUNdIrIudai9wHP+9veD39Cz8W43TpmD3BERG60ZBARWe3xt9izFW72PpDBMFuMojC8WXkMSBSRvehW5JuCbP9Z4EIR2Yl2Ma0Isv3NwPdFZAewBvhmiPL9BchDt1G3eQ9wi4jYrajtca/fQLuktgBtIX6OwRAU02bcYIhDROQG4BqlVFhHshoMM8HEKAyGOENEfoqem3xlrGUxGMAoCoMhZojIfwA3ei2+Syn16VjIYzD4w7ieDAaDwRAQE8w2GAwGQ0CMojAYDAZDQIyiMBgMBkNAjKIwGAwGQ0D+Pynrzr/CUdmWAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "clustered_data['sst'] -= clustered_data['sst'].mean(dim='anchor_year')\n", + "\n", + "clustered_data['sst'].sel(cluster_labels=1).plot.line(x='anchor_year')\n", + "clustered_data['sst'].sel(cluster_labels=-1).plot.line(x='anchor_year')\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Selecting different intervals for the precursor data (i.e. lags) will have a significant impact on the clusters that will be found:" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[array([-1., 0.]),\n", + " array([-2., -1., 0., 1.]),\n", + " array([-1., 0., 1.]),\n", + " array([-1., 0., 1.]),\n", + " array([-2., -1., 0.]),\n", + " array([-2., -1., 0.]),\n", + " array([-3., -2., -1., 0.]),\n", + " array([-1., 0.]),\n", + " array([-1., 0.]),\n", + " array([-1., 0.]),\n", + " array([0.]),\n", + " array([-1., 0.])]" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "n_intervals = field_resampled.i_interval.size\n", + "cdata = [0]*n_intervals\n", + "for interval in range(n_intervals):\n", + " cdata[interval] = s2spy.rgdr._rgdr.reduce_to_clusters(\n", + " field_resampled.sel(i_interval=interval), eps_km=600, alpha=0.05)\n", + "\n", + "[c.cluster_labels.values for c in cdata]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3.10.5 ('ai4s2s')", + "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.10.5" + }, + "orig_nbformat": 4, + "vscode": { + "interpreter": { + "hash": "a87ac8a6238ac8fffcd070769ea0ffc634a0ba702c52b706c3554dda7777b3fd" + } + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/s2spy/rgdr/__init__.py b/s2spy/rgdr/__init__.py index e8e3421..76e2a89 100644 --- a/s2spy/rgdr/__init__.py +++ b/s2spy/rgdr/__init__.py @@ -1,4 +1,3 @@ """Response Guided Dimensionality Reduction.""" -from . import _map_analysis -from . import _map_regions +from . import _rgdr from ._rgdr import RGDR diff --git a/s2spy/rgdr/_map_analysis.py b/s2spy/rgdr/_map_analysis.py deleted file mode 100644 index f20194b..0000000 --- a/s2spy/rgdr/_map_analysis.py +++ /dev/null @@ -1,88 +0,0 @@ -"""Spatial-temporal data analysis utils. - -A toolbox for spatial-temporal data analysis, including regression, -correlation, auto-correlation and relevant utilities functions. -""" -import numpy as np -import xarray as xr -from scipy.stats import pearsonr as _pearsonr - - -def _pearsonr_nan(x: np.ndarray, y: np.ndarray): - """NaN friendly implementation of scipy.stats.pearsonr. Calculates the correlation - coefficient between two arrays, as well as the p-value of this correlation. - - Args: - x: 1-D array - y: 1-D array - Returns: - r_coefficient - p_value - - """ - if np.any(np.isnan(x)) or np.any(np.isnan(y)): - return np.nan, np.nan - return _pearsonr(x, y) - - -def correlation(field: xr.DataArray, target: xr.DataArray, corr_dim: str = "time"): - """Calculate correlation maps. - - Args: - field: Spatial data with a dimension named `corr_dim`, over which each - location should have the Pearson correlation coefficient calculated with the - target data. - target: Data which has to be correlated with the spatial data. Requires a - dimension named `corr_dim`. - - Returns: - r_coefficient: DataArray filled with the correlation coefficient for each - non-`corr_dim` coordinate. - p_value: DataArray filled with the two-tailed p-values for each computed - correlation coefficient. - """ - assert ( - corr_dim in target.dims - ), f"input target does not have contain the '{corr_dim}' dimension" - assert ( - corr_dim in field.dims - ), f"input field does not have contain the '{corr_dim}' dimension" - assert np.all( - [dim in field.dims for dim in target.dims] - ), "Field and target dims do not match" - - return xr.apply_ufunc( - _pearsonr_nan, - field, - target, - input_core_dims=[[corr_dim], [corr_dim]], - vectorize=True, - output_core_dims=[[], []], - ) - - -def partial_correlation(field, target, z): - """Calculate partial correlation maps.""" - raise NotImplementedError - - -def regression(field, target): - """Regression analysis on entire maps. - - Methods include Linear, Ridge, Lasso. - """ - raise NotImplementedError - - -def save_map(): - """Save calculated coefficients. - - Store calculated coefficients and significance values, and - save them as netcdf files. - """ - raise NotImplementedError - - -def load_map(): - """Load coefficients from saved netcdf maps.""" - raise NotImplementedError diff --git a/s2spy/rgdr/_map_regions.py b/s2spy/rgdr/_map_regions.py deleted file mode 100644 index 791b101..0000000 --- a/s2spy/rgdr/_map_regions.py +++ /dev/null @@ -1,24 +0,0 @@ -"""Regions clustering utils. - -A module for clustering regions based on the given correlation -between spatial fields and target timeseries. -""" - -def spatial_mean(): - """Calculate 1-d timeseries for each precursor region. Precursor regions - are integer label masks within the np.ndarray labels. - """ - raise NotImplementedError - -def cluster_dbscan(): - """Perform DBSCAN clustering from vector array or distance matrix. - - Density-Based Spatial Clustering of Applications with Noise (DBSCAN). - Clusters gridcells together which are of the same sign and in proximity to - each other using DBSCAN. - - The input data will be processed in this function to ensure that the distance - is free of the impact from spherical curvature. The actual Euclidean distance - will be obtained and passed to the DBSCAN clustering function. - """ - raise NotImplementedError diff --git a/s2spy/rgdr/_rgdr.py b/s2spy/rgdr/_rgdr.py index 0aa8121..ff89f25 100644 --- a/s2spy/rgdr/_rgdr.py +++ b/s2spy/rgdr/_rgdr.py @@ -1,23 +1,22 @@ """Response Guided Dimensionality Reduction.""" +from typing import Tuple +import numpy as np +import xarray as xr +from scipy.stats import pearsonr as _pearsonr +from sklearn.cluster import DBSCAN class RGDR: """Response Guided Dimensionality Reduction.""" + def __init__(self, timeseries, lag_shift: int = 0): """Instantiate an RGDR operator.""" self.timeseries = timeseries self.lag_shift = lag_shift - def lag_shifting(self): - """loop through data with lag shift upto the given `lag_shift` value.""" - # To do: loop through lags - # To do: call `map_analysis` and `map_regions` and manage - # inputs and outputs. - raise NotImplementedError - def _map_analysis(self): """Perform map analysis. - + Use chosen method from `map_analysis` and perform map analysis. """ raise NotImplementedError @@ -33,3 +32,264 @@ def _clustering_regions(self): def fit(self, data): """Perform RGDR calculations with given data.""" raise NotImplementedError + + def transform(self, data): + """Apply RGDR on data, based on the fit model""" + raise NotImplementedError + + +radius_earth_km = 6371 +surface_area_earth_km2 = 5.1e8 + + +def spherical_area(latitude: float, dlat: float, dlon: float = None) -> float: + """Approximate the area of a square grid cell on a spherical (!) earth. + Returns the area in square kilometers of earth surface. + + Args: + latitude (float): Latitude at the center of the grid cell (deg) + dlat (float): Latitude grid resolution (deg) + dlon (float): Longitude grid resolution (deg), optional in case of a square grid. + + Returns: + float: Area of the grid cell (km^2) + """ + if dlon is None: + dlon = dlat + dlon = np.radians(dlon) + dlat = np.radians(dlat) + + lat = np.radians(latitude) + h = np.sin(lat + dlat / 2) - np.sin(lat - dlat / 2) + spherical_area = h * dlon / np.pi * 4 + + return spherical_area * surface_area_earth_km2 + + +def cluster_area(ds: xr.Dataset, cluster_label: float) -> float: + """Determines the total area of a cluster. Requires the input dataset to have the + variables `area` and `cluster_labels`. + + Args: + ds (xr.Dataset): Dataset containing the variables `area` and `cluster_labels`. + cluster_label (float): The label for which the area should be calculated. + + Returns: + float: Area of the cluster `cluster_label`. + """ + return ( + ds["area"].where(ds["cluster_labels"] == cluster_label).sum(skipna=True).values + ) + + +def remove_small_area_clusters(ds: xr.Dataset, min_area_km2: float) -> xr.Dataset: + """Removes the clusters where the area is under the input threshold. + + Args: + ds (xr.Dataset): Dataset containing `cluster_labels` and `area`. + min_area_km2 (float): The minimum allowed area of each cluster + + Returns: + xr.Dataset: The input dataset with the labels of the clusters set to 0 when the + area of the cluster is under the `min_area_km2` threshold. + """ + clusters = np.unique(ds["cluster_labels"]) + areas = [cluster_area(ds, c) for c in clusters] + valid_clusters = np.array([c for c, a in zip(clusters, areas) if a > min_area_km2]) + + ds["cluster_labels"] = ds["cluster_labels"].where( + np.isin(ds["cluster_labels"], valid_clusters), 0 + ) + + return ds + + +def weighted_groupby(ds: xr.Dataset, groupby: str, weight: str, method: str = "mean") -> xr.Dataset: + """Apply a weighted reduction after a groupby call. xarray does not currently support + combining `weighted` and `groupby`. An open PR adds supports for this functionality + (https://github.com/pydata/xarray/pull/5480), but this branch was never merged. + + Args: + ds (xr.Dataset): Dataset containing the coordinates or variables specified in + the `groupby` and `weight` kwargs. + groupby (str): Coordinate which should be used to make the groups. + weight (str): Variable in the Dataset containing the weights that should be used. + method (str): Method that should be used to reduce the dataset, by default + 'mean'. Supports any of xarray's builtin methods, e.g. median, min, max. + + Returns: + xr.Dataset: Dataset reduced using the `groupby` coordinate, using weights = + based on `ds[weight]`. + """ + groups = ds.groupby(groupby) + + # find stacked dim name: + group_dims = list(groups)[0][1].dims # Get ds of first group + stacked_dims = [dim for dim in group_dims.keys() if "stacked_" in dim] + + reduced_data = [ + getattr(g.weighted(g[weight]), method)(dim=stacked_dims) for _, g in groups + ] + return xr.concat(reduced_data, dim=groupby) + + +def masked_spherical_dbscan( + ds: xr.Dataset, alpha: float = 0.05, eps_km: float = 600, min_area_km2: float = None +) -> xr.Dataset: + """Determines the clusters based on sklearn's DBSCAN implementation. Alpha determines + the mask based on the minimum p_value. Grouping can be adjusted using the `eps_km` + kwarg. Cluster labels are negative for areas with a negative correlation coefficient + and positive for areas with a positive correlation coefficient. Areas without any + significant correlation are put in the cluster labelled '0'. + + Args: + ds (xr.Dataset): Dataset containing 'latitude' and 'longitude' dimensions in + degrees. Must also contain 'p_val' and 'corr' to base the groups on. + alpha (float): Value below which the correlation is significant enough to be + considered + eps_km (float): The maximum distance (in km) between two samples for one to be + considered as in the neighborhood of the other. This is not a maximum bound + on the distances of points within a cluster. This is the most important + DBSCAN parameter to choose appropriately. + min_area_km2 (float): The minimum area of a cluster. Clusters smaller than this + minimum area will be discarded. + + Returns: + xr.Dataset: Dataset grouped by the DBSCAN clusters. + """ + ds = ds.stack(coord=["latitude", "longitude"]) + coords = np.asarray(list(ds["coord"].values)) # turn array of tuples to 2d-array + coords = np.radians(coords) + + # Prepare labels, default value is 0 (not in cluster) + labels = np.zeros(len(coords)) + + for sign, sign_mask in zip([1, -1], [ds["corr"] >= 0, ds["corr"] < 0]): + mask = np.logical_and(ds["p_val"] < alpha, sign_mask) + + if np.sum(mask) > 0: # Check if the mask contains any points to cluster + db = DBSCAN( + eps=eps_km / radius_earth_km, + min_samples=1, + algorithm="auto", + metric="haversine", + ).fit(coords[mask]) + + labels[mask] = sign * (db.labels_ + 1) + + ds["cluster_labels"] = ("coord", labels) + + ds = ds.unstack(("coord")) + + dlat = np.abs(ds.latitude.values[1] - ds.latitude.values[0]) + dlon = np.abs(ds.longitude.values[1] - ds.longitude.values[0]) + ds["area"] = spherical_area(ds.latitude, dlat, dlon) + + if min_area_km2: + ds = remove_small_area_clusters(ds, min_area_km2) + return ds + + +def reduce_to_clusters( + ds: xr.Dataset, alpha: float = 0.05, eps_km: float = 600, min_area_km2: float = None +) -> xr.Dataset: + """Perform DBSCAN clustering on a prepared Dataset, and then group the data by their + determined clusters, taking the weighted mean. The weight is based on the area of + each grid cell. + + Density-Based Spatial Clustering of Applications with Noise (DBSCAN). + Clusters gridcells together which are of the same sign and in proximity to + each other using DBSCAN. + + The input data will be processed in this function to ensure that the distance + is free of the impact from spherical curvature. The actual geodesic distance + will be obtained and passed to the DBSCAN clustering function. + + Args: + ds: Dataset prepared to have p_val and corr. + alpha (float): Value below which the correlation is significant enough to be + considered + eps_km (float): The maximum distance (in km) between two samples for one to be + considered as in the neighborhood of the other. This is not a maximum bound + on the distances of points within a cluster. This is the most important + DBSCAN parameter to choose appropriately. + + """ + ds = masked_spherical_dbscan( + ds, alpha=alpha, eps_km=eps_km, min_area_km2=min_area_km2 + ) + + ds = weighted_groupby(ds, groupby="cluster_labels", weight="area") + + return ds + + +def _pearsonr_nan(x: np.ndarray, y: np.ndarray) -> Tuple[float, float]: + """NaN friendly implementation of scipy.stats.pearsonr. Calculates the correlation + coefficient between two arrays, as well as the p-value of this correlation. However, + instead of raising an error when encountering NaN values, this function will return + both the correlation coefficient and the p-value as NaN. + + Args: + x: 1-D array + y: 1-D array + Returns: + r_coefficient + p_value + + """ + if np.any(np.isnan(x)) or np.any(np.isnan(y)): + return np.nan, np.nan + return _pearsonr(x, y) + + +def correlation( + field: xr.DataArray, target: xr.DataArray, corr_dim: str = "time" +) -> Tuple[xr.DataArray, xr.DataArray]: + """Calculate correlation maps. + + Args: + field: Spatial data with a dimension named `corr_dim`, over which each + location should have the Pearson correlation coefficient calculated with the + target data. + target: Data which has to be correlated with the spatial data. Requires a + dimension named `corr_dim`. + corr_dim: Dimension over which the correlation coefficient should be calculated. + + Returns: + r_coefficient: DataArray filled with the correlation coefficient for each + non-`corr_dim` coordinate. + p_value: DataArray filled with the two-tailed p-values for each computed + correlation coefficient. + """ + assert ( + corr_dim in target.dims + ), f"input target does not have contain the '{corr_dim}' dimension" + assert ( + corr_dim in field.dims + ), f"input field does not have contain the '{corr_dim}' dimension" + assert np.all( + [dim in field.dims for dim in target.dims] + ), "Field and target dims do not match" + + return xr.apply_ufunc( + _pearsonr_nan, + field, + target, + input_core_dims=[[corr_dim], [corr_dim]], + vectorize=True, + output_core_dims=[[], []], + ) + + +def partial_correlation(field, target, z): + """Calculate partial correlation maps.""" + raise NotImplementedError + + +def regression(field, target): + """Regression analysis on entire maps. + + Methods include Linear, Ridge, Lasso. + """ + raise NotImplementedError diff --git a/setup.cfg b/setup.cfg index f932077..6db6c40 100644 --- a/setup.cfg +++ b/setup.cfg @@ -32,6 +32,7 @@ zip_safe = False include_package_data = True packages = find: install_requires = + netcdf4 numpy pandas scikit-learn diff --git a/tests/test_rgdr/test_data/sst_daily_1979-2018_5deg_Pacific_175_240E_25_50N.nc b/tests/test_rgdr/test_data/sst_daily_1979-2018_5deg_Pacific_175_240E_25_50N.nc new file mode 100644 index 0000000..9f02596 Binary files /dev/null and b/tests/test_rgdr/test_data/sst_daily_1979-2018_5deg_Pacific_175_240E_25_50N.nc differ diff --git a/tests/test_rgdr/test_data/tf5_nc5_dendo_80d77.nc b/tests/test_rgdr/test_data/tf5_nc5_dendo_80d77.nc new file mode 100644 index 0000000..36a6a2d Binary files /dev/null and b/tests/test_rgdr/test_data/tf5_nc5_dendo_80d77.nc differ diff --git a/tests/test_rgdr/test_map_analysis.py b/tests/test_rgdr/test_map_analysis.py deleted file mode 100644 index 9525dbd..0000000 --- a/tests/test_rgdr/test_map_analysis.py +++ /dev/null @@ -1,87 +0,0 @@ -"""Tests for the s2s._RGDR.map_analysis module.""" -import numpy as np -import pandas as pd -import pytest -import xarray as xr -from s2spy.rgdr import _map_analysis - - -# pylint: disable=protected-access - - -class TestMapAnalysis: - @pytest.fixture(autouse=True) - def dummy_dataarray(self): - time = pd.date_range("20161001", "20211001", freq="60d") - - da = xr.DataArray( - np.tile(np.arange(len(time)), (2, 2, 1)), - dims=["lat", "lon", "time"], - coords={ - "time": time, - "lat": np.arange(0, 2), - "lon": np.arange(0, 2), - }, - ) - return da - - @pytest.fixture(autouse=True) - def dummy_timeseries(self, dummy_dataarray): - return dummy_dataarray.isel(lat=0, lon=0).drop_vars(["lat", "lon"]) - - def test_pearsonr(self): - result = _map_analysis._pearsonr_nan([0, 0, 1], [0, 1, 1]) - expected = (0.5, 0.66666667) - np.testing.assert_array_almost_equal(result, expected) - - def test_pearsonr_nan(self): - result = _map_analysis._pearsonr_nan([np.nan, 0, 1], [0, 1, 1]) - expected = (np.nan, np.nan) - np.testing.assert_array_almost_equal(result, expected) - - def test_correlation(self, dummy_dataarray, dummy_timeseries): - c_val, p_val = _map_analysis.correlation(dummy_dataarray, dummy_timeseries) - - np.testing.assert_equal(c_val.values, 1) - np.testing.assert_equal(p_val.values, 0) - - def test_correlation_dim_name(self, dummy_dataarray, dummy_timeseries): - da = dummy_dataarray.rename({"time": "i_interval"}) - ts = dummy_timeseries.rename({"time": "i_interval"}) - c_val, p_val = _map_analysis.correlation( - da, ts, corr_dim="i_interval" - ) - - np.testing.assert_equal(c_val.values, 1) - np.testing.assert_equal(p_val.values, 0) - - def test_correlation_wrong_target_dim_name(self, dummy_dataarray, dummy_timeseries): - ts = dummy_timeseries.rename({"time": "dummy"}) - with pytest.raises(AssertionError): - _map_analysis.correlation(dummy_dataarray, ts) - - def test_correlation_wrong_field_dim_name(self, dummy_dataarray, dummy_timeseries): - da = dummy_dataarray.rename({"time": "dummy"}) - with pytest.raises(AssertionError): - _map_analysis.correlation(da, dummy_timeseries) - - def test_correlation_wrong_target_dims(self, dummy_dataarray): - ts = dummy_dataarray.rename({'lat': 'latitude'}) - with pytest.raises(AssertionError): - _map_analysis.correlation(dummy_dataarray, ts) - - def test_partial_correlation(self): - with pytest.raises(NotImplementedError): - _map_analysis.partial_correlation("field", "target", "z") - - def test_regression(self): - with pytest.raises(NotImplementedError): - _map_analysis.regression("field", "target") - - def test_save_map(self): - with pytest.raises(NotImplementedError): - _map_analysis.save_map() - - def test_load_map(self): - with pytest.raises(NotImplementedError): - _map_analysis.load_map() diff --git a/tests/test_rgdr/test_map_regions.py b/tests/test_rgdr/test_map_regions.py deleted file mode 100644 index d25252a..0000000 --- a/tests/test_rgdr/test_map_regions.py +++ /dev/null @@ -1,16 +0,0 @@ -"""Tests for the s2s._RGDR.map_regions module.""" -import pytest -from s2spy.rgdr import _map_regions - - -# pylint: disable=protected-access - - -class TestMapRegions: - def test_spatial_mean(self): - with pytest.raises(NotImplementedError): - _map_regions.spatial_mean() - - def test_cluster_dbscan(self): - with pytest.raises(NotImplementedError): - _map_regions.cluster_dbscan() diff --git a/tests/test_rgdr/test_rgdr.py b/tests/test_rgdr/test_rgdr.py index 248b4c3..bf646ba 100644 --- a/tests/test_rgdr/test_rgdr.py +++ b/tests/test_rgdr/test_rgdr.py @@ -1,6 +1,11 @@ """Tests for the s2s._RGDR.rgdr module.""" +import numpy as np +import pandas as pd import pytest +import xarray as xr from s2spy import RGDR +from s2spy.rgdr import _rgdr +from s2spy.time import AdventCalendar # pylint: disable=protected-access @@ -8,21 +13,152 @@ class TestRGDR: """Test RGDR methods.""" + @pytest.fixture(autouse=True) def dummy_rgdr(self): - rgdr = RGDR('timeseries') - return rgdr + return RGDR("timeseries") def test_init(self): - rgdr = RGDR('timeseries') + rgdr = RGDR("timeseries") assert isinstance(rgdr, RGDR) - def test_lag_shifting(self, dummy_rgdr): - with pytest.raises(NotImplementedError): - dummy_rgdr.lag_shifting() - def test_fit(self, dummy_rgdr): with pytest.raises(NotImplementedError): dummy_rgdr.fit("data") - + def test_transform(self, dummy_rgdr): + with pytest.raises(NotImplementedError): + dummy_rgdr.transform("data") + + +class TestCorrelation: + @pytest.fixture(autouse=True) + def dummy_dataarray(self): + time = pd.date_range("20161001", "20211001", freq="60d") + + return xr.DataArray( + np.tile(np.arange(len(time)), (2, 2, 1)), + dims=["lat", "lon", "time"], + coords={ + "time": time, + "lat": np.arange(0, 2), + "lon": np.arange(0, 2), + }, + ) + + @pytest.fixture(autouse=True) + def dummy_timeseries(self, dummy_dataarray): + return dummy_dataarray.isel(lat=0, lon=0).drop_vars(["lat", "lon"]) + + def test_pearsonr(self): + result = _rgdr._pearsonr_nan([0, 0, 1], [0, 1, 1]) + expected = (0.5, 0.66666667) + np.testing.assert_array_almost_equal(result, expected) + + def test_pearsonr_nan(self): + result = _rgdr._pearsonr_nan([np.nan, 0, 1], [0, 1, 1]) + expected = (np.nan, np.nan) + np.testing.assert_array_almost_equal(result, expected) + + def test_correlation(self, dummy_dataarray, dummy_timeseries): + c_val, p_val = _rgdr.correlation(dummy_dataarray, dummy_timeseries) + + np.testing.assert_equal(c_val.values, 1) + np.testing.assert_equal(p_val.values, 0) + + def test_correlation_dim_name(self, dummy_dataarray, dummy_timeseries): + da = dummy_dataarray.rename({"time": "i_interval"}) + ts = dummy_timeseries.rename({"time": "i_interval"}) + c_val, p_val = _rgdr.correlation(da, ts, corr_dim="i_interval") + + np.testing.assert_equal(c_val.values, 1) + np.testing.assert_equal(p_val.values, 0) + + def test_correlation_wrong_target_dim_name(self, dummy_dataarray, dummy_timeseries): + ts = dummy_timeseries.rename({"time": "dummy"}) + with pytest.raises(AssertionError): + _rgdr.correlation(dummy_dataarray, ts) + + def test_correlation_wrong_field_dim_name(self, dummy_dataarray, dummy_timeseries): + da = dummy_dataarray.rename({"time": "dummy"}) + with pytest.raises(AssertionError): + _rgdr.correlation(da, dummy_timeseries) + + def test_correlation_wrong_target_dims(self, dummy_dataarray): + ts = dummy_dataarray.rename({"lat": "latitude"}) + with pytest.raises(AssertionError): + _rgdr.correlation(dummy_dataarray, ts) + + def test_partial_correlation(self): + with pytest.raises(NotImplementedError): + _rgdr.partial_correlation("field", "target", "z") + + def test_regression(self): + with pytest.raises(NotImplementedError): + _rgdr.regression("field", "target") + + +test_file_path = "./tests/test_rgdr/test_data" + + +class TestMapRegions: + @pytest.fixture(autouse=True) + def dummy_calendar(self): + return AdventCalendar(anchor_date=(8, 31), freq="30d") + + @pytest.fixture(autouse=True) + def example_target(self): + return xr.open_dataset(f"{test_file_path}/tf5_nc5_dendo_80d77.nc").sel( + cluster=3 + ) + + @pytest.fixture(autouse=True) + def example_field(self): + return xr.open_dataset( + f"{test_file_path}/sst_daily_1979-2018_5deg_Pacific_175_240E_25_50N.nc" + ) + + @pytest.fixture(autouse=True) + def example_corr(self, dummy_calendar, example_target, example_field): + field = dummy_calendar.resample(example_field) + target = dummy_calendar.resample(example_target) + + field["corr"], field["p_val"] = _rgdr.correlation( + field.sst, target.ts.sel(i_interval=0), corr_dim="anchor_year" + ) + + return field + + @pytest.fixture(autouse=True) + def expected_labels(self): + return np.array( + [ + [0.0, 0.0, 0.0, 0.0, -1.0, -1.0, 0.0, -2.0, -2.0, -2.0, 0.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -2.0, -2.0, -2.0, 0.0, 0.0], + [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -2.0, -2.0, -2.0, 0.0, 0.0], + [0.0, 0.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -2.0, 0.0, 0.0], + [1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, -2.0, -2.0, -2.0, -2.0, -2.0], + ] + ) + + def test_dbscan(self, example_corr, expected_labels): + clusters = _rgdr.masked_spherical_dbscan(example_corr.sel(i_interval=1)) + + np.testing.assert_array_equal(clusters["cluster_labels"], expected_labels) + + def test_dbscan_min_area(self, example_corr, expected_labels): + clusters = _rgdr.masked_spherical_dbscan( + example_corr.sel(i_interval=1), min_area_km2=3000**2 + ) + expected_labels[expected_labels == -1] = 0 # Small -1 cluster is missing + + np.testing.assert_array_equal(clusters["cluster_labels"], expected_labels) + + def test_cluster(self, example_corr): + clustered_data = _rgdr.reduce_to_clusters( + example_corr.sel(i_interval=1), eps_km=600, alpha=0.04 + ) + + np.testing.assert_array_equal(clustered_data["cluster_labels"], [-1, 0, 1]) + + assert set(clustered_data.sst.dims) == {"anchor_year", "cluster_labels"}