diff --git a/docs/source/notebooks/paper/Get_ERA5_data_and_Calibrate_a_hydrological_model.ipynb b/docs/source/notebooks/paper/Get_ERA5_data_and_Calibrate_a_hydrological_model.ipynb new file mode 100644 index 00000000..732c43ad --- /dev/null +++ b/docs/source/notebooks/paper/Get_ERA5_data_and_Calibrate_a_hydrological_model.ipynb @@ -0,0 +1,695 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Create a watershed model" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Hydrological models typically need geographical information about watersheds being simulated: latitude and longitude, area, mean altitude, land-use, etc. This notebook shows how to obtain this information using remote services that are made available for users in PAVICS-Hydro. These services connect to a digital elevation model (DEM) and a land-use data set to extract relevant information.\n", + "\n", + "The DEM used in the following is the [EarthEnv-DEM90](https://www.earthenv.org/DEM), while the land-use dataset is the [North American Land Change Monitoring System](http://www.cec.org/tools-and-resources/north-american-environmental-atlas/north-american-land-change-monitoring-system). Other data sources could be used, given their availability through the Web Coverage Service (WCS) protocol.\n", + "\n", + "Since these computations happen on a specific Geoserver hosted in PAVICS, we need to establish a connection to that service. While the steps are a bit more complex, the good news is that you only need to change a few items in this notebook to taylor results to your needs. For example, this first code snippet is boilerplate and should not be changed.\n", + "\n", + "We will also setup a hydrological model, calibrate it, and save the parameters for future use.\n", + "\n", + "We will be using the Mistassini river as the test-case for this example, but you can substitute the data for any catchment of your liking. We provide:\n", + "\n", + "1- Streamflow observations (Water Survey Canada station 02RD003)\n", + "\n", + "2- Watershed boundaries in the form of shapefiles (all shape files .shp, .shx, .prj, .dbf, etc. zipped into a single file. The platform will detect and unzip the file to extract the required data)\n", + "\n", + "The rest will be done by PAVICS-Hydro." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "
\n", + "\n", + "\n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + " \n", + "
featuresNameOfficialIDFlagPAVICSSourceAreageometry
01MISTASSINI (RIVIERE) EN AMONT DE LA RIVIERE MI...02RD0031HYDAT9870POLYGON ((-72.26250 48.87917, -72.27720 48.881...
\n", + "
" + ], + "text/plain": [ + " features Name OfficialID \\\n", + "0 1 MISTASSINI (RIVIERE) EN AMONT DE LA RIVIERE MI... 02RD003 \n", + "\n", + " FlagPAVICS Source Area geometry \n", + "0 1 HYDAT 9870 POLYGON ((-72.26250 48.87917, -72.27720 48.881... " + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 1, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAJcAAAD5CAYAAAAqYadTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAeTklEQVR4nO2deXhU5d33P7+Z7CQhCYHsJAHDvhOQRVwBFVvFrUpF8cGntFZal6eLVq/3tdZWq7UufdzR6utWrXWhoiJVURDZ952QsGSDLJCVJCRzv3/MBLJMkpnMnDkzZ87nunLNzJn73OeXyTf3uedefl9RSmFiogUWvQMwMS6muEw0wxSXiWaY4jLRDFNcJpphistEM0JcKSQih4AaoAVoVkrlisj1wIPAcGCyUmpjF+deBjwNWIElSqlHe7peYmKiysrKciU0E53ZtGlTuVKqv7P3XBKXg4uUUuVtXu8ErgFe7OoEEbECzwKzgEJgg4gsVUrt7u5CWVlZbNzoVKsmfoaIHO7qPXfE1Q6l1B5H5d0VmwzkKaXyHWX/AVwFdCsuE2Pgap9LAV+IyCYRWeRG/WnA0TavCx3HOiEii0Rko4hsLCsrc+MSJv6Kq+KarpSaAFwO3CEi57t4nrNmzel8k1LqJaVUrlIqt39/p7dwkwDDJXEppYodj8eBD7Hf7lyhEMho8zodKHYnQJPApUdxiUgfEYlpfQ7Mxt6Zd4UNQI6IZItIGHAjsLS3wZoEFq60XEnAahHZBqwHlimlPheRq0WkEJgKLBOR5QAikioinwIopZqBxcByYA/wnlJqlxa/iIn/If645CY3N1eZQxGBgYhsUkrlOnvPHKE30Yxej3OZGJsdhVU8sWIfIRYLsZEhzJ+SyYSB8W7VYYoriGlsbqHFpogKay+Dt9cd4f8u3cnplrNdpjCrxRSXietU1DZxweNfM2FgPDlJ0dQ1tnCyvomv97UfxB6ZGst9lw93u35TXEFMTEQIp1sU6woqWVdQ2WW5vOO1fJ9fwWWjkt2q3+zQBzGRoVaXyjU22/jTp3twd2TBFFcQ8589x1wu+8OxKT0tUuiEKa4gRSnF89/ku1R2SFI0C6ZmuX0Ns88VhNhsivc3F7Lt6Mkey144tD9P3zCevlGhbl/HFFcQcbiijvUFlXy+s5Qv9x7vtmxGQiTP/ngCY9Ljen09U1xBxE1L1lF44pRLZS8ZluSRsMAUV1BR39TSY5lrxqcxNiOO3Cz3BkydYYoriEiNi6CyrqnL96+dkM5j143BanHvW2FXmN8Wg4hhybHdvn/luFSvCQtMcRmaFpui6tRpWmyKusZm8o7Xdlv+p29s5EhFvdeub94WA5CPtxYxJj2O7MQ+Z44ppfhwSxGnW2zcMGkgAH/5Yh/PrzwIQFiIhaZmW7f1Npy2sWR1Pg9dNcorcZri8lOqG06zYtcx9h+rYc3BCgbEhPPYdWPYePgEd7+7lfdvn0Y2dnGtL6jk2a/z+GZ/GRaB+KgwbIozwgJ6FFYrB8u6b93cwRSXH6KU4pFP9/DO+qPtjl/7/BpKqhqwKVi5r4zhybG8vf4If/jk7DZQm4JfvLOFUKv+PR5TXH6EUoqPtxbz4rf57Cmp7vT+oTb9oWe+PMBr3xVQ3dDcqVxjs41GF1uqjrTYvLfs3RSXH7HqQDl3vbvV5fLOhOUpM4cnea0uU1x+wKbDJ/hsRwmf7SzVOxQuH53itbpMcfkBr685xNJt+u8VnpyVQFpcpNfq07/XZ0JEqP5/hrioUP56w1iv1mm2XDrT2NzC0UrXJpO9TXp8JAunZ6OAcRl9SY+P8mr9prh0Zn1BJd/nV+hy7bnj0lh4XrZm9evfHgc50wcnMn/KQF2u/enOEo5VN2hWvykunbFYhOsmZvRcUAPyy+p4f1OhZvWb4vIDPtdxCGLDoa63lHmKKS6d2Xb0JC+vcm2jhBZsOnwCmxdH5dtiduh1oqnZxgvfHOR/v8rz6pSLu9Q1NlNZ30RidLjX6zbFpQM2m+Lu97aybHuJ3qFw7YR0TYQFprg0p6ymkbKaRmIiQoiNCGXzkRM89Z/9bCus0js0wkMs3DN7iGb1m+LSkPLaRuY8s4qymka9Q3FKmNVCRIhrW/p7g9mh15DE6HB+MkO7QUpPqWlsZuPhE+2OKaVYue8485esY3vhSY/qN1suDTndYmNHUed1Wf7E3pJqRqbGEhcVyqurC3hr3RFKqho4Z0A0CX3CPKrbFJdGNLfYWPz2Zpbvcj3Zhx48sWI/T6zYT5jVQlPL2QWGd16S4/FcoykuDThUXscv3tnCjiL9O+2u0lZYACFe2GJm9rk0YFdxdUAJqyNpcZGMyYjzuB6z5dKAmIjA/Vijwqz882dTSfXCokFP/BYTgHeBLOAQ8COl1AlXzvU4aj/FZlNsL6rizbVdusT5PVeOTfXaoKpLJgcOgeS29VsUkceASqXUoyJyLxCvlPqtK+f2RCCaHKw6UMav/rmNY9X+OablDhaBlL6RDEyI4m8/Ht+t2LQyObgKeN3x/HVgrgd1BTwTM+O5alwaXky1oBs2BUUnT/F9fgV/XbG/1/V44reYpJQqAXA8DnDjXMMRFRbC7+YMJzo8cPtbznh73RHye7kL29VPYrpSqlhEBgArRGSvG9fodK5S6tuOhRzCWwQwcKA+KzM95VRTiyZ7CfVkYmZ8u5wU7uCJ3+IxEUkBcDw6zYPoqlejEcw8tVwyrBe/vnSo21mcW/HEb3EpsMBRbAHwsRvnGobSqrOCKjWYuMZlxHFudkKvz++13yLwKDBLRA4Asxyv2/ktdnOuIThYVsvdju33tY3NPPLpHn0D8jKzRiT1utUC02/RIx79bC8vfHOQYckxlNc2Ul7bdUrIQCTMamHfw5d1KzDTb9EFPthcSEF5ncvlj1c38JZjsHRvaY3hhAUQHmrxqOUKWnEppaiobWR9QSXPrzzI//xzG9WnTrt8/kOf7Kam0VjfDDvS22+JrRhrUMZFGk638F9/39Bpp3Odi2I5UdfEJ36w/l1rfnFxjkfnB424dhZV8c3+MmIiQvh2f7nTLfQn6s+2XDabYum2Yi4blUxEB3cvZ4nZjER8VChXjElhdFpfj+oJCnHtKq5i3strqelhgPOe97by5d5jRIZa2XzkJHtKqnkz+lzOy0lsV263wcUVHRHCw3NHe1yP4cW1u7iaG19c61L/qLHZxgebi9ode/yLfUEnrp7+CV3FsOI6WlnP7z7cwdr8inZeze5S7mTnzu5i44tLKeXRN0UwsLieW5nHqgMur/LpEovF/s3y1OkW3lx7mHX5lew/VuOFCP2XFpuivqmFPh5OwhtSXGvyyjul2e4tRytP8eSK/XyyvYR8N8bBAp1/bS7kll4YeLbFcOJSSvHnz91ZtNEzz3yV59X6AoHjXlj0aLhB1H9tLvKLrfKBztTB/Tyuw3ji0jCZWTBxygVvxp4wnLjunOnZqLKJnRYvLGgwnLgmZyWQ1c+7WYmDkdxMz51iDScui0WYPyVT7zACmklZ8fTzwvYyw4kL4IZJGUzzQoc0GJk7LpV3fjLFK3UZUlwxEaH8/MJz9A4joAixCLedl81Dc0cR4iU7PcONcwFU1DbyyGfGWnKsNS8vyOWioV3tDuwdhmu56puaWfjaBnYZfP7Pm4jAsOQYr9drGHEppXj523yueW6NOYjqJkrBj178nuKT3vUgMoy4lm4r5o+f7mFvqbEnlbXiaOUpXvzmYM8F3cAQ4mqxKe791w69wwh4IsO82wU3hLisFmFEaqzeYQQ8Ni9vMzTMt8WZw5PYdLhTejATFxid1peZw5O4aJh30ygYRlzDU7z/bSdYePz6MQxL9n7Lb4jbIni/SQ8mkmMjNKnXMOKKiQjVO4SAZFhyDHFRnuWb7wrDiMukd2jpmGYYcS3R0bMwkMkrq6WmwfU0Bu5gCHEppSjy8uhysHBO/2jNsiEaQlwiwuKLzBWo7nLbedksv+t80ryQc94ZhhAXQHLfCK9YigQTk7ISsGj4mRlGXOMy4vj6VxfqHUZAcde7W/hqr3bGV4YZRAWIjTSHI1wlMtTK/CkDGZ6i3bSZocTljyk4/ZEQi/DZnTPI8jC5W08Y5rYIEBcVppkZuJG4cmyq5sICg4kLtFlRaTSG+OgzMpy4+kaZ/a7usFqEwf2jfXItw4nLW4nLjMoPx6Qwa0SST65lOHHVajSVYRRyknzXbXBJXCJySER2iMhWEdnoOJYgIitE5IDj0en+bxG5TET2iUiew5dRU567aSI/mZGt9WUCFl/dEsG9lusipdS4Nm4J9wJfKqVygC8dr9shIlbgWeByYAQwT0RGeBhztyT3jeC+y4ebo/UdiI0I4c/XjubSkb65JYL2Zp6TgTylVL5Sqgn4h+M8TbFYhAEx5pBEWx6/fiw3TBrocZ5Td9DazDMNaJs/stBxrBMiskhENorIxrKyMhfD6pqkvtqsrgxU/rhsDw2nPc+55Q6uimu6UmoC9tvbHSJyvovnOfs3cTqM7m2/xUGJvutbBAJHKuvZUeTbzcJam3kWAhltXqcDxZ4E7Cq/vnQofc25xnb4OsW5pmaewAYgR0SyRSQMuNFxnuYk943gx+cGpp2xVpwzwLetuSsT10nAh46OYAjwtlLqcxHZALwnIrcBR4DrwW7mCSxRSs1RSjWLyGJgOWAFXlVK7dLiF+nIR1uKvL49PdBparb59Ho9iksplQ+MdXK8ArjEyfFiYE6b158Cn3YspzXLd5Wi4d6DgOS+D3aw/K7zfTZFZrgR+lZCvZTAzEiUVjdwot53pqOG/AuUVJ3iyz3arbAMVAb17+OTpTatGFJcD/17N3VeyKNuNHJ83KE3nLhO1DWx4VCl3mH4JZ6ac7qL4cQV3yeM1b+9mO0PzjYzOrchPT6S/54xyKfXNJy4ACJCrcRGhPLna8cQZnbsASg8cYqfvrHJZR9vb2DoTz4jIYoLhno351Qg883+Mt7b6B2rQFcwtLgAoj00pDQSEaEWny5FMvwnb67qOsu/bp/GyFTfdeoN33JplWQjEPl8Z6lPr2d4cc3ISdQ7BL9hzcEKn17P8OK6ZkIaoVbz5gjwuznDfXo9w4srJiKUqYPN1gvwuQ+l4cUF8KerR3Hf5cN8vp7J36j38ZRYUIgrPT6Kn14wmAkD4/QORVciw6w+vV5QiAugqv40X+11thI7ePjN+9t9umAwaMT1xIp9lNf6bi2TP/LV3uP85Yt9Prte0Ijr1mlZ9OujTb71QOKjLUXYfLREN2jENah/NO/fPk3vMHTneE0j2320xSxoxAWYu7CxW7Gkx2uTvbkjQSWu8BALf7p6NNk+XOrrT4RYhGdvmuCz7ItBJa4Qq4V5kzP44ZgUvUPRhdkjk5iY6TQZkSYYflVEKyv32b8pVdQ2UVLVoHc4uvCfPcd59us8Fp0/yCe7o4Km5fpoSxE7i6qDVlhg3xT76uoCqk75JkGeIcXl7Ku2t/2bA5U/zB3lsz6X4T7xfaU1LHpjI9MGJ3L1+DTqm5qpbmim8ES93qHpSqhVmD0imctGJvvsmoYS1+GKOha8up7S6gYOVxzhnfVH9A7Jb3j8urHMHe80NZpmGEZcRyvruf6F7zle06h3KH7HzVMyfS4sMEifq6ymkRtfWmsKqwsuGKLPDqiAb7mOVzdw7wc7TDPPbvBhGtR2BHzLta6gMuiX0vTE7W9tZm2+b9fPQ4CLq7axmceX+24JSaDS1Gzjxy+v5eZX1vHhlkLqm3yzIyqgb4uvrznEkcrgHmJwFZuCVQfKWXWgnKiwnbx480Rm5GjbFwtoca0vMLPZ9ES/PmFU1NkXSQ5LjmFESiwZCVE+yaERsOJafaBcl35EoBAeYmH+lEx+NXsoe0urqaxr4uJhA3xqchAQ4lJK8chne5k6uB8XDbV7Kfzli300+jiBbCBx8bAB3DNrCJFhVsYP9N1KiLYERId+zcEKXvo2nz8u24PNplBKcbCsVu+w/JrPdpYy+8lv2VPi29zzbQkIcSVGh9MnzEre8Vru/2gnH24pMn0VXaDo5Cn+/l2BbtcPCHENTY7hyRvGAfDO+iPc8942fQMKID7bUepzz59WXBaXiFhFZIuIfOJ4PVZEvnf4MP5bRGK7OK+TV2NvmD0ymV/NHtLb04OGjku4axqbWblPn0Fmd1quO4E9bV4vAe5VSo3G7gf0627O7ejV2CvuuOgcJmXp0zkNBMKsFpYuns4bt00+c2xSVjxTBumTG9ZVp9h04ArsgmplKPCt4/kK4FrvhuY0DiZlJWh9mYClqcXGM18eYEZOfxZMzeT2Cwfz+sLJxEXps1/T1aGIp4DfAG0NkncCV2I3lLqe9u5kbWn1alTAi0qpl5wVcvg4LgIYOLBrQ6ixGXEuhhycvLK6gCvHpvH7q0bpHYpLrmU/AI4rpTZ1eGshdu/FTdhF19VeeZe8Gl31W7xk2ACGJAV3tprusCm494PtNLfoPwboym1xOnCliBzCbiN8sYi8qZTaq5SarZSaCLwDOLUI68KrsdeEWC1MM/Ntdcuu4mpW7vPcbddTehSXUuo+pVS6UioLu1/iV0qp+SIyAEBELMADwAsdz+3Gq9EjDlXUeVqF4fGHPQOejHPNE5H9wF7s7q9/B7vfooi0WuAlAatFZBuwHlimlPrck4CVUhw4Zo7O90SxH2yhc2tuUSm1EljpeP408LSTMmf8FrvyavSE3SXV5qpTF/gur5xD5XU+dSnrSECM0LcloU8YP79wMCNTY/Fhvv6AY1dxNb94Z4vP0iU5I+DEldI3kt9cNoxlv5zB9gcvNU3Su2FHURWvmnOLvSM6PMSniTUCkbd13LsZ0OICTHH1QH5ZnW7fHANeXFeMTqG/mdStW55ccUCX6wa8uLIS+/DB7dNIi/NNtrxA5ItdpWw67Pv9BgEvLrD7Ko5J963FbiBR09jMtc9/z77SGp9e1xDiamq28b25WaNHnvnyADuLqlDKzObsMsdrGnRbbRlILNtRwg/+tpr7P/J4Bs4lDCGu9Pgo/jZvgjmo2g19wqyk9I2gT5jVTP7mLrNGJLHh/pnsKalhd0kVu4ur+WR7Cc06jlD7E9POSeTlWzxaCOw2hmi5WukXHc55OYksOn8wPxiTagqrDUUnfD8fayhxtfL9wQp+9mbHtY3Bze6SagrKfbtUyZDi2lVcZbZaTvC1x7Vh+lxtMTfMnuVXs4cwJCmGZpvyuXOIIcUVE2HIX8ttYiNC+O8Zg4gI9a2JZyuGvC3OGpFEiDkuwdzxaboJCwzacmX268MDVwznwX/v1jsUXfj8rhlEhFiJ19lf0pAtF8D8KZkMTYrpuaABSY2LJCuxj+4LKQ0rrhCrhSULcrlk2IAzxwYl9iEm3JCNdTvK/CRluqE/6YyEKJ6bP4FPd5QwObsfaXGR7CmpZsGr6w2ds375rlJ+fuE5eodh3JarlfAQK1ePTz+z3mt4Siz3zDJ2tpwPNxf5bOVDdxheXM447Qdb3bXkwPFa9pT4du2WM4JSXMl9jbdqNS4qlFCrMHN4EtmJffj2gP7b+Q3d5+qKzH5ReofgVVL7RrDmvktQSiEiKKUor+0qL4zvCMqWKyM+in46jwF5k+Ep9qSOrWnARcQvNq0Epbgiw6y8cdu5xEUZY0Ntjp+O5wWluABGpMby7qKphmjB3l53mHV+uIcgaMUF7bNEBzLVDc28sfaw3mF0IqjFBfhF38QbxET43y0+6MV1yiC7hj7ZVszWoyf1DqMdQS+u5hb9R7K9QU1jMze/so6qU6f1DuUMQS+uiFDjfAQ1Dc28sipf7zDOYJxPtpdY9DKA1oiXVxX4fCNGVwS9uIzm2XjqdAsz//oNd/1jC3nH9Z1fDGpxKaX436/z9A7D67TYFB9tLebq59awZFU+tY36bFgJanGJCKeajPFt0Rk1Dc08vGwPlzyxksIT9TQ2+/Z3DWpxAcQGQU7VY9WNzF+yjoWvbfBplsGgFtf+YzVnWq7YiBCGp8QyKs39LNGPXjOaBVMzdV+z3h2HKur5Lq+CW15Zz1vrDlPqgzz1Li+5ERErsBEoUkr9QETGYnfNiAYOATcppTp53orIZdjz1VuBJUqpR70RuDd4eNmeM/2R/5k9lAXTsgCoqG3kq73HeXlVPvtdMFSobWzm91eN4qYpmVz97HfUeXCrvWTYANYVVGrWT8ovr+P+D3dyPztJj49kUlYCEzPjmZSVQM6AaCxe3JKnqd+iQ5DPYjeVGoHddWNE78P1LjkDzhpUtS5bAXtCk+tzM1i6+DxuOrdrB7VW9joy9g1JiuHpG8fT29GNYckxPDNvPP9YNMUnG3sLT5ziwy1FPPDRTi596lsWv7OZD7cUsq+0xivLpLX2W5wM5Cml8pVSTdiNqa7qfbjeZd7kgWd2Bw1L6bxsJSLUysNzRzEjp3sjq4vb7DCaOSKJhdOzexXP3tIaahubGZXWl9f+axLRTnYqabl76dMdpdz97jYufepbjlR63jdzteV6CrvfYtvF561+i9C132IacLTN60LHsU6IyCIR2SgiG8vKfLNE95wB0bxy6yT+vfg8YruY+BURnrxhHAOcTHBHh4eQFhfJhUPbW/gtOn8QoVZ78xUZauXyUcm8MH9ij32y8BAL/R2J2SZmJvDSzRMZmRrLRUP7c+OkDO6eOYT198/kozumc+u0LE2XC204dMLjOnr8N2jrtygiF7Z5ayHwjIj8H2Apzv0Wnd0gnLa3DpPPlwByc3N9OuE3uodkvYnR4dx2XjaPfLa33fEnbxjHrBFJnconxUbw/E0T6R8TzojUWEKtFpbvKu1x3i89PrJdn2faOYks++WMTuXGZcQxLiOO+68Yzr7SGg6W1XKwrI78slpWHSj3yvzipsOVXDcx3aM6XGljW/0W5wARQKzDb3E+dos7RGQI9ttmRwpp36KlY3c4CzhumZrFK6sL2u137K7vO7ON6N5ed4QHPtrR4zWi3Vw2E2q1MCqtL6PSzv5z3LRkLd/leT7r4I2WS1O/RWADkCMi2SIS5jh/qcdR60BkmJWrx9vv6FaLMC4jjoEJrm30KKtppKt0YWlxkVw8bABzRidzbrbn/t25mWfr8GTadGRqbM+FesCT3uE8EbnD8fwD2vgtYh9ymKOUahaRxcBy7EMRryqldnkUsY7MmzyQqYP7kZuV4LSz3RU/OT+b51bm0dhsY2RqLE3NNlpsilunZzFv8kBCrd4bbrxrZg4Lp2cTHRFCZV0Tn2wv5vOdpWw4VNmlwJ1xwZCuraBdRfxhZ25HcnNz1caNG/UOw6u89l0B4wfG62YAv7OoilteXU9lnWtbzgb178Obt51Lag/OJCKySSnlNJNvUI/Q+5Jbp2frJiyAUWl9+fzOGTxwxXDCQyw8du2YdmN7Hckvq+Pa59fw9+8KqG/q3YCu2XIFIWU1jfSPCedwRR0XPL6yx/JxUaF8cPs0BvWP7vSe2XKZtKN1U0pmvz4MTYphREosv79yZJdfUE7Wn+a+D3a47ToblNv5Tc7y5A3jGBAbTmJ0OPOnZPLs13n8dcX+TuXWFVRysKzWrQ24priCnBFthhysFuGXl+Qwb/JASqsaeGvdYT7aWkTDafvETHltEzmdx4y7xBSXSSf6x4TTPyacR9PHcN+c4fy/NYdYvrvU7Y692aE38QizQ2+iC6a4TDTDFJeJZpjiMtEMU1wmmmGKy0QzTHGZaIYpLhPNMMVlohl+OUIvImWAXkk+E4Fyna7tDXwdf6ZSyumyVb8Ul56IyMaupjMCAX+K37wtmmiGKS4TzTDF1ZmX9A7AQ/wmfrPPZaIZZstlohmmuEw0IyjFJSLvishWx88hEdnqOD65zfFtInJ1F+cniMgKETngeIz3k/hnicgmEdnheLy4i/MfFJGiNnXM0STOYO9zicgTQJVS6iERiQKaHGkIUoBtQKpSqrnDOY8BlUqpR0XkXiBeKfVb30ffKf7xwDGlVLGIjAKWK6U6pawSkQeBWqXUX7SMLShbrlbE7n75I+AdAKVUfRshRdBFuifsCexedzx/HZirYZhd4iT+LUqp1ixCu4AIEdHNOSuoxQXMwP6ffqD1gIicKyK7gB3Azzq2Wg6SlFIlAI7HAU7K+IJO8bfhWmCLUqrRyXsAi0Vku4i8qtVt3bC3RRH5D5Ds5K37lVIfO8o8jz2t5hNOzh+OvVU6XynV0OG9k0qpuDavTyilvPoH8iR+ERmJPVXVbKXUQSd1J2Gff1TAH4AUpdRCb8YP2F0kgvEH+57NY0B6N2W+BnKdHN/n+IMApAD7/CV+7An29gPTXawnC9ipRYzBfFucCexVShW2HnAkqQtxPM/EnlT4kJNzlwILHM8XAB9rG6pTnMUfBywD7lNKfdfViY4vK61cjT2/rffRuwXR6wd4DXufqu2xm7F3hLcCm4G5bd5bgqMVA/oBXwIHHI8JfhL/A0CdI/7WnwFO4n8De59yO/Z/lBQtYjRsn8tEf4L5tmiiMaa4TDTDFJeJZpjiMtEMU1wmmmGKy0QzTHGZaMb/B+s1iCOAzdTGAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# We need to import a few packages required to do the work\n", + "import json\n", + "import os\n", + "import rasterio\n", + "import warnings\n", + "import geopandas as gpd\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import rioxarray as rio\n", + "import xarray as xr\n", + "from birdy import WPSClient\n", + "from zipfile import ZipFile\n", + "\n", + "# The platform provides lots of user warnings and information points. We will disable them for now.\n", + "warnings.filterwarnings('ignore')\n", + "\n", + "# This is the URL of the Geoserver that will perform the computations for us.\n", + "url = os.environ.get(\"WPS_URL\", \"https://pavics.ouranos.ca/twitcher/ows/proxy/raven/wps\")\n", + "\n", + "# Connect to the PAVICS-Hydro Raven WPS server\n", + "wps = WPSClient(url)\n", + "\n", + "# Name of the watershed boundaries file that is uploaded to the server\n", + "feature_url = 'shapefile_basin_574_HYSETS.zip'\n", + "\n", + "# Prepare a plot of the catchment\n", + "df = gpd.read_file(feature_url)\n", + "display(df)\n", + "df.plot()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Generic watershed properties\n", + "\n", + "Now that we have delineated a watershed, lets find the zonal statistics and other properties using the `shape_properties` process. This process requires a `shape` argument defining the watershed contour, the exterior polygon.\n", + "\n", + "Once the process has completed, we extract the data from the response, as follows. Note that you do not need to change anything here. The code will work and return the desired results." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'id': '0',\n", + " 'features': 1,\n", + " 'Name': 'MISTASSINI (RIVIERE) EN AMONT DE LA RIVIERE MISTASSIBI',\n", + " 'OfficialID': '02RD003',\n", + " 'FlagPAVICS': 1,\n", + " 'Source': 'HYDAT',\n", + " 'Area': 9870,\n", + " 'area': 9569366229.461872,\n", + " 'centroid': [-72.7431067594341, 49.848278236356585],\n", + " 'perimeter': 727186.8545423769,\n", + " 'gravelius': 2.0970051622223544}" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "{'area': 9569.366229461872,\n", + " 'longitude': -72.7431067594341,\n", + " 'latitude': 49.848278236356585,\n", + " 'gravelius': 2.0970051622223544,\n", + " 'perimeter': 727186.8545423769}" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "shape_resp = wps.shape_properties(shape=feature_url)\n", + "\n", + "[properties,] = shape_resp.get(asobj=True)\n", + "prop = properties[0]\n", + "display(prop)\n", + "\n", + "area = prop[\"area\"] / 1000000.0\n", + "longitude = prop[\"centroid\"][0]\n", + "latitude = prop[\"centroid\"][1]\n", + "gravelius = prop[\"gravelius\"]\n", + "perimeter = prop[\"perimeter\"]\n", + "\n", + "shape_info = {\n", + " \"area\": area,\n", + " \"longitude\": longitude,\n", + " \"latitude\": latitude,\n", + " \"gravelius\": gravelius,\n", + " \"perimeter\": perimeter,\n", + "}\n", + "display(shape_info)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note that these properties are a mix of the properties of the original file where the shape is stored, and properties computed by the process (area, centroid, perimeter and gravelius). Note also that the computed area is in m², while the \"SUB_AREA\" property is in km², and that there are slight differences between the two values due to the precision of HydroSHEDS and the delineation algorithm.\n", + "\n", + "## Land-use information\n", + "\n", + "Now we extract the land-use properties of the watershed using the `nalcms_zonal_stats` process. As mentionned, it uses the [North American Land Change Monitoring System](http://www.cec.org/tools-and-resources/north-american-environmental-atlas/north-american-land-change-monitoring-system) dataset, and retrieve properties over the given region.\n", + "\n", + "With the `nalcms_zonal_stats_raster` process, we also return the grid with variable accessors (`gdal`, `rasterio`, or `rioxarray`) depending on what libraries are available in our runtime environment (The following examples show `rioxarray`-like access)." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "stats_resp = wps.nalcms_zonal_stats_raster(\n", + " shape=feature_url, \n", + " select_all_touching=True, \n", + " band=1, \n", + " simple_categories=True\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Here we will get the raster data and compute statistics on it. It is also possible to download the extracted raseter offline (please see the tutorial for the steps on how to do this)." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Downloading to /tmp/tmp6bmgq3k0/subset_1.tiff.\n" + ] + }, + { + "data": { + "text/plain": [ + "'Land use ratios'" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "{'Ocean': 0.0,\n", + " 'Forest': 0.7246596208414477,\n", + " 'Shrubs': 0.14616312094792794,\n", + " 'Grass': 0.04322426804857576,\n", + " 'Wetland': 0.013300924493021603,\n", + " 'Crops': 0.00395034960218003,\n", + " 'Urban': 0.0035571063310866975,\n", + " 'Water': 0.06514460973576021,\n", + " 'SnowIce': 0.0}" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "'Land use percentages'" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "text/plain": [ + "{'Ocean': '0.0 %',\n", + " 'Forest': '72.47 %',\n", + " 'Shrubs': '14.62 %',\n", + " 'Grass': '4.32 %',\n", + " 'Wetland': '1.33 %',\n", + " 'Crops': '0.4 %',\n", + " 'Urban': '0.36 %',\n", + " 'Water': '6.51 %',\n", + " 'SnowIce': '0.0 %'}" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "features, statistics, grid0 = stats_resp.get(asobj=True)\n", + "lu = statistics[0]\n", + "total = sum(lu.values())\n", + "\n", + "land_use = {k: (v / total) for (k, v) in lu.items()}\n", + "display(\"Land use ratios\", land_use)\n", + "\n", + "land_use_pct = {\n", + " k: \"{} %\".format(np.round(v/total*100, 2))\n", + " for (k, v) in lu.items()\n", + "}\n", + "display(\"Land use percentages\", land_use_pct)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Terrain information from the DEM\n", + "\n", + "Here we collect terrain data, such as elevation, slope and aspect, from the DEM. We will do this using the `terrain_analysis` WPS service, which by default uses DEM data from [EarthEnv-DEM90](https://www.earthenv.org/DEM).\n", + "\n", + "Note here that while the feature outline is defined above in terms of geographic coordinates (latitude, longitude), the DEM is projected onto a 2D cartesian coordinate system (here NAD83, the Canada Atlas Lambert projection). This is necessary to perform slope calculations. For more information on this, see: https://en.wikipedia.org/wiki/Map_projection\n", + "\n", + "The DEM data returned in the process response here shows `rioxarray`-like access but using the URLs we can open the files however we like." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "{'elevation': 423.6657935442332,\n", + " 'slope': 3.949426174669343,\n", + " 'aspect': 148.5591531205915}" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "terrain_resp = wps.terrain_analysis(\n", + " shape=feature_url, \n", + " select_all_touching=True, \n", + " projected_crs=3978\n", + ")\n", + "\n", + "properties, dem0 = terrain_resp.get(asobj=True)\n", + "\n", + "elevation = properties[0][\"elevation\"]\n", + "slope = properties[0][\"slope\"]\n", + "aspect = properties[0][\"aspect\"]\n", + "\n", + "terrain = {\n", + " \"elevation\": elevation, \n", + " \"slope\": slope, \n", + " \"aspect\": aspect\n", + "}\n", + "display(terrain)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Overview\n", + "\n", + "A synthesis of all watershed properties can be created by merging the various dictionaries created. This allows users to easily access any of these values, and to provide them to a Raven model as needed." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "pycharm": { + "name": "#%%\n" + } + }, + "outputs": [ + { + "data": { + "text/plain": [ + "{'area': 9569.366229461872,\n", + " 'longitude': -72.7431067594341,\n", + " 'latitude': 49.848278236356585,\n", + " 'gravelius': 2.0970051622223544,\n", + " 'perimeter': 727186.8545423769,\n", + " 'Ocean': 0.0,\n", + " 'Forest': 0.7246596208414477,\n", + " 'Shrubs': 0.14616312094792794,\n", + " 'Grass': 0.04322426804857576,\n", + " 'Wetland': 0.013300924493021603,\n", + " 'Crops': 0.00395034960218003,\n", + " 'Urban': 0.0035571063310866975,\n", + " 'Water': 0.06514460973576021,\n", + " 'SnowIce': 0.0,\n", + " 'elevation': 423.6657935442332,\n", + " 'slope': 3.949426174669343,\n", + " 'aspect': 148.5591531205915}" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "all_properties = {**shape_info, **land_use, **terrain}\n", + "display(all_properties)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Getting input data\n", + "\n", + "Now that we have all the geographic information for our watershed, we can get the input meteorological data required to calibrate and run the model. Here we use an in-house solution that keeps updated ERA5 reanalysis datasets available with little to no wait." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "import datetime as dt\n", + "import s3fs\n", + "import fsspec\n", + "import intake\n", + "from clisops.core import subset\n", + "\n", + "\n", + "# This will be our input section, where we control what we want to extract. \n", + "# We know which watershed interests us, it is the shapefile we used to get geographic data:\n", + "\n", + "basin_contour = 'shapefile_basin_574_HYSETS.zip' # Contains .shp, .shx and other shapefile extension files\n", + "\n", + "# Also, we can specify which timeframe we want to extract. Here let's focus on a 10-year period\n", + "reference_start_day = dt.datetime(1980, 12, 31)\n", + "reference_stop_day = dt.datetime(1991, 1, 1) \n", + "# Notice we are using one day before and one day after the desired period of 1981-01-01 to 1990-12-31.\n", + "# This is to account for any UTC shifts that might require getting data in a previous or later time.\n", + "\n", + "\n", + "# Get the ERA5 data from the Wasabi/Amazon S3 server. \n", + "catalog_name = 'https://raw.githubusercontent.com/hydrocloudservices/catalogs/main/catalogs/atmosphere.yaml'\n", + "cat=intake.open_catalog(catalog_name)\n", + "ds=cat.era5_reanalysis_single_levels.to_dask()\n", + "\n", + "'''\n", + "Get the ERA5 data. We will rechunk it to a single chunck to make it compatible with other codes on the platform, especially bias-correction.\n", + "We are also taking the daily min and max temperatures as well as the daily total precipitation.\n", + "'''\n", + "# We will add a wrapper to ensure that the following operations will preserve the original data attributes, such as units and variable names.\n", + "with xr.set_options(keep_attrs=True):\n", + " \n", + " ERA5_reference=subset.subset_shape(ds.sel(time=slice(reference_start_day,reference_stop_day)), basin_contour)\n", + " ERA5_tmin=ERA5_reference['t2m'].resample(time='1D').min().chunk(-1,-1,-1)\n", + " ERA5_tmax=ERA5_reference['t2m'].resample(time='1D').max().chunk(-1,-1,-1)\n", + " ERA5_pr=ERA5_reference['tp'].resample(time='1D').sum().chunk(-1,-1,-1)\n", + "\n", + " # Change the units\n", + " ERA5_tmin = ERA5_tmin - 273.15 # K to °C\n", + " ERA5_tmin.attrs['units'] = 'degC'\n", + "\n", + " ERA5_tmax = ERA5_tmax - 273.15 # K to °C\n", + " ERA5_tmax.attrs['units'] = 'degC'\n", + "\n", + " ERA5_pr = ERA5_pr * 1000 # m to mm\n", + " ERA5_pr.attrs['units'] = 'mm'\n", + " \n", + " # Average the variables spatially\n", + " ERA5_tmin = ERA5_tmin.mean({'latitude','longitude'})\n", + " ERA5_tmax = ERA5_tmax.mean({'latitude','longitude'})\n", + " ERA5_pr = ERA5_pr.mean({'latitude','longitude'})\n", + " \n", + " # Ensure that the precipitation is non-negative, which can happen with some reanalysis models.\n", + " ERA5_pr[ERA5_pr < 0] = 0\n", + " \n", + " # Transform them to a dataset such that they can be written with attributes to netcdf\n", + " ERA5_tmin = ERA5_tmin.to_dataset(name='tmin',promote_attrs=True)\n", + " ERA5_tmax = ERA5_tmax.to_dataset(name='tmax',promote_attrs=True)\n", + " ERA5_pr = ERA5_pr.to_dataset(name='pr',promote_attrs=True) \n", + "\n", + " # Write to disk.\n", + " ERA5_tmin.to_netcdf('ERA5_tmin.nc')\n", + " ERA5_tmax.to_netcdf('ERA5_tmax.nc')\n", + " ERA5_pr.to_netcdf('ERA5_pr.nc')\n", + " " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Set-up the hydrological model\n", + "\n", + "Now that we have geographic and meteorological input data available, we can setup a Raven hydrological model and calibrate it. Many more details can be found in the documentation and tutorial notebooks. Note that the cell will take a bit of time to run as the model is being calibrated. " + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Nash-Sutcliffe value is: [0.6901]\n", + "GR4JCN.Params(GR4J_X1=0.05489502, GR4J_X2=5.117167, GR4J_X3=555.5967, GR4J_X4=4.135658, CEMANEIGE_X1=25.37618, CEMANEIGE_X2=0.8970489)\n", + "[5.489502e-02 5.117167e+00 5.555967e+02 4.135658e+00 2.537618e+01\n", + " 8.970489e-01]\n" + ] + } + ], + "source": [ + "# Import the template for our model, which will be the GR4J model with the CemaNeige snow model in our case. \n", + "from ravenpy.models import GR4JCN, GR4JCN_OST\n", + "\n", + "# Combine the input and streamflow data into a single element that will be passed to Raven\n", + "forcing = ('ERA5_tmax.nc', 'ERA5_tmin.nc', 'ERA5_pr.nc', 'Qobs_574_HYSETS.nc')\n", + "\n", + "# Using Ostrich with the GR4JCN model. Start by creating the calibration model\n", + "model = GR4JCN_OST()\n", + "\n", + "# Create the HRU for the watershed\n", + "hru = GR4JCN.LandHRU(area = all_properties['area'],\n", + " elevation = all_properties['elevation'],\n", + " latitude = all_properties['latitude'],\n", + " longitude = all_properties['longitude'])\n", + "\n", + "# Establish the start date for the calibration\n", + "import datetime as dt\n", + "start_date=dt.datetime(1981, 1, 1)\n", + "end_date = dt.datetime(1985,12,31)\n", + "\n", + "# Starting point parameters. We can provide this to help accelerate the calibration process, but they are not strictly necessary.\n", + "params = (0.529, -3.396, 407.29, 1.072, 16.9, 0.053)\n", + "\n", + "#lower and upper bounds for the parameters. Note that there are 6 values, each corresponding to the GR4JCN parameter in that position.\n", + "lower = (0.01, -15.0, 10.0, 0.0, 1.0, 0.0)\n", + "upper = (2.5, 10.0, 700.0, 7.0, 30.0, 1.0)\n", + "\n", + "# Optimization algorithm. Multiple options are available, see OSTRICH documentation for more information. Here, DDS is used as it is powerful and\n", + "# particularly useful for optimizations with small evaluation budgets. See: \n", + "#\n", + "# Tolson, B.A. and Shoemaker, C.A., 2007. Dynamically dimensioned search algorithm for computationally efficient watershed model calibration. Water \n", + "# Resources Research, 43(1)\n", + "#\n", + "# for more details.\n", + "algorithm = 'DDS' \n", + "\n", + "# Maximum number of model evaluations. We only use 200 here to keep the computation time as low as possible, but you will want to increase this \n", + "# for operational use, perhaps to 500 or 1000.\n", + "max_iterations = 200\n", + "\n", + "# Random seed. We will provide one for consistency purposes, but operationnaly this should not be provided.\n", + "random_seed=0\n", + "\n", + "# Here is where we launch the model calibration with the desired parameters and options. The process builds the model in the background and will return optimized parameters.\n", + "model(\n", + " ts = forcing,\n", + " hrus = (hru,),\n", + " start_date=start_date,\n", + " end_date=end_date,\n", + " params=params,\n", + " lowerBounds=lower,\n", + " upperBounds=upper,\n", + " algorithm=algorithm,\n", + " random_seed=random_seed, # Remove this for operational use!\n", + " max_iterations=max_iterations,\n", + " overwrite=True,\n", + ")\n", + "\n", + "# Get the model diagnostics including Nash-Sutcliffe Efficiency\n", + "d = model.diagnostics\n", + "\n", + "# Print the NSE and the parameter set in 2 different ways:\n", + "print('Nash-Sutcliffe value is: ' + str(d['DIAG_NASH_SUTCLIFFE']))\n", + "print(model.calibrated_params) # With explanations of what these parameters are\n", + "print(model.optimized_parameters) # Just the array that could be used in another process. This is what people will typically want to use." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Run the calibrated hydrological model on a validation period\n", + "\n", + "Now that the hydrological model has been calibrated, we can use these parameters to run the model on an independent period for validation" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "defaultdict(list,\n", + " {'observed data series': ['HYDROGRAPH_ALL'],\n", + " 'filename': ['Qobs_574_HYSETS.nc'],\n", + " 'DIAG_NASH_SUTCLIFFE': [0.485837],\n", + " 'DIAG_RMSE': [136.853]})" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "# Instantiate a GR4JCN model\n", + "model_validation = GR4JCN()\n", + "\n", + "# Run the model with validation period dates and data, as well as calibrated parameters. Notice that the forcing data stays the same: Raven will select the data of the correct dates from the forcing data.\n", + "model_validation(\n", + " ts=forcing,\n", + " start_date=dt.datetime(1986,1,1),\n", + " end_date=dt.datetime(1990,12,31),\n", + " hrus=(hru,), # Careful how this must be passed! This is due to the capability of running in distributed mode as well.\n", + " params=model.optimized_parameters,\n", + " overwrite=True, # OPTIONAL: We can do this to overwrite old files with the new ones generated in this run (output files, etc.)\n", + ")\n", + "\n", + "# plot the diagnostics so we can see the validation Nash-Sutcliffe Efficiency:\n", + "model_validation.diagnostics" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Save parameters and properties for future use\n", + "\n", + "We can now save the calibrated parameters and basin properties for use in the second example, where we will look at the impacts of climate change." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "np.savetxt('optimized_parameters_GR4JCN_basin_574_HYSETS.txt',model.optimized_parameters)\n", + "json.dump(all_properties, open(\"properties_basin_574_HYSETS.txt\",'w'))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.8.13" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/docs/source/notebooks/paper/Perform_a_climate_change_impact_assessment.ipynb b/docs/source/notebooks/paper/Perform_a_climate_change_impact_assessment.ipynb new file mode 100644 index 00000000..f7b06a9b --- /dev/null +++ b/docs/source/notebooks/paper/Perform_a_climate_change_impact_assessment.ipynb @@ -0,0 +1,484 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "12ea4c8a-ac6a-4ad0-b0d7-1aa3dbf39fa9", + "metadata": {}, + "source": [ + "# Evaluate the impacts of climate change on streamflow\n", + "\n", + "This notebook shows how we can work with climate model data to evaluate the impacts of climate change on hydrology. We will apply a processing chain to the data from the previous example on the Mistassini catchment and we will see how to implement the full \n", + "process in PAVICS-Hydro.\n", + "\n", + "We will first start by setting up the required information to make a GR4JCN model, including calibrated parameters and properties generated in the previous example: " + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "01c1784a-5527-4a9f-9829-d8593aecf204", + "metadata": {}, + "outputs": [], + "source": [ + "# Import required packages\n", + "import xclim\n", + "import gcsfs\n", + "import intake\n", + "import json\n", + "import warnings\n", + "import numpy as np\n", + "import xarray as xr\n", + "import datetime as dt\n", + "import xclim.sdba as sdba\n", + "import datetime as dt\n", + "from ravenpy.models import GR4JCN\n", + "from clisops.core import subset, average\n", + "\n", + "# The platform provides lots of user warnings and information points. We will disable them for now.\n", + "warnings.filterwarnings('ignore')\n", + "\n", + "# Load the properties and calibrated parameters from the previous example\n", + "all_properties = json.load(open(\"properties_basin_574_HYSETS.txt\"))\n", + "parameters = np.loadtxt('optimized_parameters_GR4JCN_basin_574_HYSETS.txt')\n", + "\n", + "# Get the basin contour shapefile, as the previous example\n", + "basin_contour = 'shapefile_basin_574_HYSETS.zip'\n", + "\n", + "# Setup the GR4JCN HRU that will be used to build the model later, making use of the previously-obtained catchment properties\n", + "hru = GR4JCN.LandHRU(area = all_properties['area'],\n", + " elevation = all_properties['elevation'],\n", + " latitude = all_properties['latitude'],\n", + " longitude = all_properties['longitude'])" + ] + }, + { + "cell_type": "markdown", + "id": "4be57177-ad3d-4fab-9f65-c7289dca5719", + "metadata": {}, + "source": [ + "## Extract climate model data\n", + "\n", + "PAVICS-Hydro makes use of the PanGEO catalog to obtain CMIP6 project data directly from Google Cloud Services. This makes obtaining data much more efficient. This next section shows how to implement this data extraction for a climate model, but you can access more climate models and datasets by following the tutorial on PAVICS-Hydro." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "98b92acb-46bb-4bfa-914c-23be01f01e3b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "historical tasmin\n" + ] + }, + { + "name": "stderr", + "output_type": "stream", + "text": [ + "hwloc/linux: Ignoring PCI device with non-16bit domain.\n", + "Pass --enable-32bits-pci-domain to configure to support such devices\n", + "(warning: it would break the library ABI, don't enable unless really needed).\n" + ] + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "historical tasmax\n", + "historical pr\n", + "ssp585 tasmin\n", + "ssp585 tasmax\n", + "ssp585 pr\n" + ] + } + ], + "source": [ + "# Start by setting up the reference and future start and end days. Reference data is used to calibrate the bias-correction algorithms that will be applied to the future data.\n", + "reference_start_day = dt.datetime(1980, 1, 1)-dt.timedelta(days=1)\n", + "reference_end_day = dt.datetime(1991, 12, 31)+dt.timedelta(days=1) # Notice we are using one day before and one day after the desired period of 1981-01-01 to 1990-12-31. This is to account for any UTC shifts that might require getting data in a previous or later time.\n", + "future_start_day = dt.datetime(2080, 1, 1)-dt.timedelta(days=1)\n", + "future_end_day = dt.datetime(2091, 12, 31)+dt.timedelta(days=1) # Notice we are using one day before and one day after the desired period of 2080-01-01 to 2090-12-31. This is to account for any UTC shifts that might require getting data in a previous or later time.\n", + "\n", + "# Climate model to use\n", + "climate_model = 'MIROC6'\n", + "\n", + "\n", + "#Get the catalog info from the pangeo dataset, which basically is a list of links to the various products.\n", + "fsCMIP = gcsfs.GCSFileSystem(token='anon', access='read_only')\n", + "col = intake.open_esm_datastore(\"https://storage.googleapis.com/cmip6/pangeo-cmip6.json\")\n", + "\n", + "# We will add a wrapper to ensure that the following operations will preserve the original data attributes, such as units and variable names.\n", + "with xr.set_options(keep_attrs=True):\n", + "\n", + " # Load the files from the PanGEO catalogs, for reference and future variables of temperature and precipitation.\n", + " out = {}\n", + " for exp in ['historical','ssp585']:\n", + " \n", + " if exp=='historical':\n", + " period_start = reference_start_day\n", + " period_end = reference_end_day\n", + " else:\n", + " period_start = future_start_day\n", + " period_end = future_end_day\n", + " \n", + " out[exp] = {}\n", + " for variable in ['tasmin','tasmax','pr']:\n", + " print(exp, variable)\n", + " query = dict(experiment_id=exp, table_id='day', variable_id=variable, member_id='r1i1p1f1', source_id=climate_model)\n", + " col_subset = col.search(require_all_on=['source_id'], **query)\n", + " mapper = fsCMIP.get_mapper(col_subset.df.zstore[0])\n", + " \n", + " # special case for precipitation, which does not have the \"height\" variable that we need to discard as for tasmax and tasmin.\n", + " if variable == 'pr':\n", + " out[exp][variable] = average.average_shape(xr.open_zarr(mapper, consolidated=True).sel(time=slice(period_start,period_end))[variable], basin_contour).chunk(-1)\n", + " else:\n", + " out[exp][variable] = average.average_shape(xr.open_zarr(mapper, consolidated=True).sel(time=slice(period_start,period_end)).reset_coords('height',drop=True)[variable], basin_contour).chunk(-1)\n", + "\n", + "# We can now extract the variables that we will need later:\n", + "historical_tasmax = out['historical']['tasmax']\n", + "historical_tasmin = out['historical']['tasmin']\n", + "historical_pr = out['historical']['pr']\n", + "future_tasmax = out['ssp585']['tasmax']\n", + "future_tasmin = out['ssp585']['tasmin']\n", + "future_pr = out['ssp585']['pr']" + ] + }, + { + "cell_type": "markdown", + "id": "68461b69-bca1-4b77-bf3e-dab11c457f23", + "metadata": {}, + "source": [ + "## Get the historical period data (observations or equivalent)\n", + "\n", + "These data were computed in the previous example, but we will re-extract them to keep the code easier to understand and integrate with other workflows." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "b08b0a67-216e-4799-9945-3a48c84809d8", + "metadata": {}, + "outputs": [], + "source": [ + "# Regenerate the ERA5 data to be sure. In an operational context, you could combine everything onto one notebook and ensure that the \n", + "# dates and locations are constant!\n", + "\n", + "# Get the ERA5 data from the Wasabi/Amazon S3 server. \n", + "catalog_name = 'https://raw.githubusercontent.com/hydrocloudservices/catalogs/main/catalogs/atmosphere.yaml'\n", + "cat=intake.open_catalog(catalog_name)\n", + "ds=cat.era5_reanalysis_single_levels.to_dask()\n", + "\n", + "'''\n", + "Get the ERA5 data. We will rechunk it to a single chunck to make it compatible with other codes on the platform, especially bias-correction.\n", + "We are also taking the daily min and max temperatures as well as the daily total precipitation.\n", + "'''\n", + "# We will add a wrapper to ensure that the following operations will preserve the original data attributes, such as units and variable names.\n", + "with xr.set_options(keep_attrs=True):\n", + " \n", + " ERA5_reference=subset.subset_shape(ds.sel(time=slice(reference_start_day,reference_end_day)), basin_contour).mean({'latitude','longitude'})\n", + " ERA5_tmin=ERA5_reference['t2m'].resample(time='1D').min().chunk(-1,-1,-1)\n", + " ERA5_tmax=ERA5_reference['t2m'].resample(time='1D').max().chunk(-1,-1,-1)\n", + " ERA5_pr=ERA5_reference['tp'].resample(time='1D').sum().chunk(-1,-1,-1)" + ] + }, + { + "cell_type": "markdown", + "id": "3e0d6625-4f71-4c97-ad2e-9672cf94693b", + "metadata": {}, + "source": [ + "## Change units\n", + "\n", + "Climate models and reanalysis datasets have often differing units to those expected by Raven. Here we update units to make them compatible" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "eb683cb1-94d1-44a9-9777-c0226cfb20f5", + "metadata": {}, + "outputs": [], + "source": [ + "# Here we need to make sure that our units are all in the correct format. You can play around with the tools we've seen thus far to explore the units\n", + "# and make sure everything is consistent.\n", + "\n", + "# Let's start with precipitation:\n", + "ERA5_pr = xclim.core.units.convert_units_to(ERA5_pr,'mm',context='hydro')\n", + "# The CMIP data is a rate rather than an absolute value, so let's get the absolute values:\n", + "historical_pr = xclim.core.units.rate2amount(historical_pr)\n", + "future_pr = xclim.core.units.rate2amount(future_pr)\n", + "\n", + "# Now we can actually convert units in absolute terms.\n", + "historical_pr = xclim.core.units.convert_units_to(historical_pr,'mm',context='hydro')\n", + "future_pr = xclim.core.units.convert_units_to(future_pr,'mm',context='hydro')\n", + "\n", + "# Now let's do temperature:\n", + "ERA5_tmin = xclim.core.units.convert_units_to(ERA5_tmin,'degC')\n", + "ERA5_tmax = xclim.core.units.convert_units_to(ERA5_tmax,'degC')\n", + "historical_tasmin = xclim.core.units.convert_units_to(historical_tasmin,'degC')\n", + "historical_tasmax = xclim.core.units.convert_units_to(historical_tasmax,'degC')\n", + "future_tasmin = xclim.core.units.convert_units_to(future_tasmin,'degC')\n", + "future_tasmax = xclim.core.units.convert_units_to(future_tasmax,'degC')" + ] + }, + { + "cell_type": "markdown", + "id": "08fafb98-77c5-4c66-ac1a-c22a39ff65f2", + "metadata": {}, + "source": [ + "## Apply bias-correction to the climate model data\n", + "\n", + "Here is where we perform the bias-correction to the reference and future climate data in order to remove biases as seen between the reference and historical data. The future dataset is then corrected with the same adjustment factors as those in the reference period." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "dc552048-2ae2-4806-aab4-344697dd1b4d", + "metadata": {}, + "outputs": [], + "source": [ + "# Use xclim utilities (sbda) to give information on the type of window used for the bias correction.\n", + "group_month_window = sdba.utils.Grouper(\"time.dayofyear\", window=15)\n", + "\n", + "# This is an adjusting function. It builds the tool that will perform the corrections.\n", + "Adjustment = sdba.DetrendedQuantileMapping.train(ref=ERA5_pr, hist=historical_pr,\n", + " nquantiles=50, \n", + " kind=\"+\", \n", + " group=group_month_window\n", + ")\n", + "\n", + "# Apply the correction factors on the reference period\n", + "corrected_ref_precip = Adjustment.adjust(historical_pr, interp=\"linear\")\n", + "\n", + "# Apply the correction factors on the future period\n", + "corrected_fut_precip = Adjustment.adjust(future_pr, interp=\"linear\")\n", + "\n", + "# Ensure that the precipitation is non-negative, which can happen with some climate models\n", + "corrected_ref_precip=corrected_ref_precip.where(corrected_ref_precip>0, 0)\n", + "corrected_fut_precip=corrected_fut_precip.where(corrected_fut_precip>0, 0)\n", + "\n", + "# Train the model to find the correction factors for the maximum temperature (tasmax) data\n", + "Adjustment = sdba.DetrendedQuantileMapping.train(ref=ERA5_tmax, hist=historical_tasmax,\n", + " nquantiles=50, \n", + " kind=\"+\", \n", + " group=group_month_window\n", + ")\n", + "\n", + "# Apply the correction factors on the reference period\n", + "corrected_ref_tasmax = Adjustment.adjust(historical_tasmax, interp=\"linear\")\n", + "\n", + "# Apply the correction factors on the future period\n", + "corrected_fut_tasmax = Adjustment.adjust(future_tasmax, interp=\"linear\")\n", + "\n", + "# Train the model to find the correction factors for the minimum temperature (tasmin) data\n", + "Adjustment = sdba.DetrendedQuantileMapping.train(ref=ERA5_tmin, hist=historical_tasmin,\n", + " nquantiles=50, \n", + " kind=\"+\", \n", + " group=group_month_window\n", + ")\n", + "\n", + "# Apply the correction factors on the reference period\n", + "corrected_ref_tasmin = Adjustment.adjust(historical_tasmin, interp=\"linear\")\n", + "\n", + "# Apply the correction factors on the future period\n", + "corrected_fut_tasmin = Adjustment.adjust(future_tasmin, interp=\"linear\")" + ] + }, + { + "cell_type": "markdown", + "id": "e90ff1e8-b473-42e6-a692-5fc4e335f7cc", + "metadata": {}, + "source": [ + "## Generate the NetCDF files\n", + "\n", + "Now that the datasets are created, we can generate files so that Raven can access them. This might take a bit of time since everything up until now has been done in a \"lazy\" framework by Python. Data processing is actually just now really starting." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "ecb69825-d452-4bf7-b8a8-34fdafe6d82c", + "metadata": {}, + "outputs": [], + "source": [ + "# Convert the reference corrected data into netCDF file. We will then apply a special code to remove a dimension in the dataset to make it applicable to the RAVEN models.\n", + "ref_dataset = xr.merge([\n", + " corrected_ref_precip.to_dataset(name=\"pr\"),\n", + " corrected_ref_tasmax.to_dataset(name=\"tasmax\"), \n", + " corrected_ref_tasmin.to_dataset(name=\"tasmin\")\n", + "])\n", + "ref_dataset.to_netcdf(\"reference_dataset_tmp.nc\")\n", + "\n", + "# Convert the future corrected data into netCDF file\n", + "fut_dataset = xr.merge([\n", + " corrected_fut_precip.to_dataset(name=\"pr\"),\n", + " corrected_fut_tasmax.to_dataset(name=\"tasmax\"), \n", + " corrected_fut_tasmin.to_dataset(name=\"tasmin\")\n", + "])\n", + "fut_dataset.to_netcdf(\"future_dataset_tmp.nc\")\n", + "\n", + "ref_dataset = xr.open_dataset(\"reference_dataset_tmp.nc\")\n", + "ref_dataset.isel(geom=0).squeeze().to_netcdf(\"reference_dataset.nc\")\n", + "\n", + "fut_dataset = xr.open_dataset(\"future_dataset_tmp.nc\")\n", + "fut_dataset.isel(geom=0).squeeze().to_netcdf(\"future_dataset.nc\")\n" + ] + }, + { + "cell_type": "markdown", + "id": "ef2fe5bb-f610-49ee-a58b-979ac7ae6459", + "metadata": {}, + "source": [ + "## Run the hydrological model with reference and future data\n", + "\n", + "Now that we have all the data we need for the reference and future periods, we can now run GR4JCN to obtain streamflow representative of each period." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "c0881a46-e0c1-4cbf-95ba-b761658b5a09", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEiCAYAAAA1YZ/LAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABnoElEQVR4nO2dd3gc1dW437OrVZcsF7lbbhiwMcbGBhNC75CEknwhENJDSOFLr+SXQBpfekggQIAACYFASEgCCb33ZoMBGzA2tjEuuBe5SJa05/fHzOzOrmZ3Z6W9q13pvs+jR7szs3PvnXLPPeWeK6qKxWKxWCyFJNLXFbBYLBZL/8MKF4vFYrEUHCtcLBaLxVJwrHCxWCwWS8GxwsVisVgsBccKF4vFYrEUHCtcLBaLxVJwrHCxGEdEVojIcQU+5zkicl8hzzlQEZFPiMgTfV2PTIjIBBFREdkhIueF/M0PRWSn+7sK03W0dMcKF0tZoqo3qeoJfV2PgYyIVInIdSKyXUTeEZGvpe2fKSLzRWSX+3+mb990EblXRDaKSNiZ3E2qerX7+0oR+Yc7cFEROcp/oKpeBOzXqwZaeoUVLhaLpaf8AJgCjAeOBr4lIieB0/kDtwM3AoOBPwO3u9sBOoBbgU/3ovwngI8A7/TiHBZDWOFiKRYHicirIrJFRK4XkWoAERksIv8VkQ3uvv+KyFjvR67JZpmItIrIchE5x7f9Cd9xKiKfE5El7nkuFxHJViHvHCLyK/c3y0XkZN/+T4rIa27Zy0Tks759R4nIKhH5loisF5G1InK6iJwiIm+IyGYR+a7v+IiIfEdE3hSRTSJyq4gMKcylTZSxl4g8KiLbXI3gb+52z6xU4Tv2ERE5N/Xncpn729dF5NgQRX4M+LGqblHV14BrgE+4+44CKoDfqmq7ql4KCHAMgKouVtVrgUU9aauq7lHV36rqE0BXT85hMYsVLpZicQ5wIjAZ2Bv4nrs9AlyPM/ptAXYDvwcQkTrgUuBkVW0ADgUWZCnjvcBBwAHAmW55uZgLLAaGAb8ArvUJpfXuORuBTwKXiMiBvt+OBKqBMcCFOJ3rR4DZwOHAhSIyyT32S8DpwJHAaGALcHlQhUSkRUS2Zvn7cIa2/Bi4D0dTGAtcFqL9/uuwzL0OFwH/zCb8RGSw246XfJtfImmK2g94WVOTF76MNVUNGKxwsRSL36vq26q6GbgYOBtAVTep6m2quktVW919R/p+Fwemi0iNqq5V1Wwj3Z+p6lZVXQk8DMwMUa+3VPUaVe3CMd2MAka4dbtTVd9Uh0dxOu7Dfb/tAC5W1Q7gFpyO+Xeq2urWcxEwwz32s8D/U9VVqtqOY1L6nyBns6quVNWmLH9/zdCWDhwhPVpV29xRfVjW42gZHar6NxyB+54sx9e7/7f5tm0DGnz7t5GKf7+ln2OFi6VYvO37/BbOqBcRqRWRq0TkLRHZDjwGNIlIVFV3Ah8CPgesFZE7RWTfLGX4be+7SHaA2Uj8RlV3uR/r3bqdLCLPuCaurcApOALEY5MrlMDRuADW+fbv9tVhPPAvT/sAXsMx54wIUcewfAvH9PSciCwSkU/l8dvVaVpG4h5lYIf7v9G3rRFo9e1vJBX/fks/xwoXS7EY5/vcAqxxP38d2AeYq6qNwBHudgFQ1XtV9XgcjeJ1HNOTcUSkCrgN+BUwQlWbgLu8evWAt3HMe34NpFpVVweU3SJO2G2mv3OCClDVd1T1M6o6GkdTukJE9gJ2uofU+g4fmfbzMWk+Kv89CiprC7AWxwTpcQBJH8oiYEbaOWfQQx+LpfywwsVSLM4XkbGuHf+7wN/c7Q04I/yt7r6LvB+IyAgROdX1vbTjjIaL5bytBKqADUCn6+jvTejzH4CLRWQ8gIg0i8hpQQe6ZrH6LH83Bf1ORD7oC4bYAijQpaobgNXAR0Qk6mo0k9N+Phz4kojEROSDwFQcYZqNG4DviROUsS/wGeBP7r5HcO7Vl8QJWf5fd/tDbl1FnKCOSvd7tSvQQ+Oet9r9Wumeo6fC31JgrHCxFIu/4vgslrl/P3G3/xaoATYCzwD3+H4TwdFs1gCbcXwxXyhGZV3/z5dwwmW3AB8G7ujFKX/n/v4+EWnFaevc3tYzjYOAZ0Vkh1vWl1V1ubvvM8A3gU04TvWn0n77LE5Y8UYcv9f/qOqmHOVdBLyJY0J7FPilqt4DTjQXTgDDx4CtwKeA093t4JgJd5PUZHbj+HnyYbH7uzHAve7n8Xmew2IIsStRWiyWUsbV9hYDbcA3VTWnaVRELgK+hqN91vl8Y5YiYYWLxWKxWApOyZjFxEkjsV5EFvq2/U1EFrh/K0Rkgbt9gojs9u37g+83s0XkFRFZKiKXWhvswEZE/pDBKf6H3L+2ZAkqODz3ry0DmZLRXETkCByH7Q2qOj1g/6+Bbar6IxGZAPw3w3HPAV/GsWnfBVyqqncbrbzFYrFYUigZzUVVH8Nx2nbD1T7OBG7Odg4RGQU0qurTbsz+DThORYvFYrEUkXJJRX04sE5Vl/i2TRSRF4HtwPdU9XGcqJFVvmNWuduyMmzYMJ0wYUIBq2uxWCz9n/nz529U1eagfeUiXM4mVWtZC7So6iYRmQ38W0T2I3iCW6DdT5x1Ic4DaGlpYd68eQWussVisfRvROStTPtKxiyWCTf30vtJTrrDzbK6yf08HyfWfm8cTWWs7+djyTDLWFWvVtU5qjqnuTlQ8FosFoulh5S8cAGOA15X1YS5y53dHHU/T8KZ/LVMVdcCrSJyiOun+RjOmhIWi8ViKSIlI1xE5GbgaWAfcdbJ8BYROovujvwjgJdF5CXgH8Dn3Gy7AJ8H/ggsxdFobKSYxWKxFJmSCUXuS+bMmaPW52KxWCz5ISLzVXVO0L6S0VwsFovF0n+wwsVisVgsBccKF4vFYrEUHCtcLBbLgGHh6m3sd+E9rG9t6+uq9HuscLFYLAOG659cwc49XTyyeENfV6XfY4WLxWIZMETdHs9GyZrHCheLxTJgiLgrcMStbDGOFS4Wi2XA4K3uFLeai3GscLFYLAMGb+1AK1vMY4WLxWIZMESs5lI0rHCxWCwDhoTPxTpdjGOFi8ViGTB4Cz5Z0WIeK1wsFsuAwfO5WMxjhYvFYhlwWJeLeaxwsVgsAwZPcbGyxTxWuFgslgGD4IUiW/FiGitcLBbLgMG6XIqHFS4Wi2XAYRUX81jhYrFYBgzJUGQrXUxjhYvFYhkwJBz6VrYYxwoXi8UyYEjkFuvjegwESka4iMh1IrJeRBb6tv1ARFaLyAL37xTfvgtEZKmILBaRE33bZ4vIK+6+S8XOmrJYLC4Js5iVLsYpGeEC/Ak4KWD7Jao60/27C0BEpgFnAfu5v7lCRKLu8VcC5wFT3L+gc1osloGIHWoWjZIRLqr6GLA55OGnAbeoaruqLgeWAgeLyCigUVWfVieQ/QbgdCMVtlgsZYt16JunZIRLFv5XRF52zWaD3W1jgLd9x6xyt41xP6dv74aInCci80Rk3oYNdj1ti2UgELHruRSNUhcuVwKTgZnAWuDX7vYg5VazbO++UfVqVZ2jqnOam5sLUFWLxVLqJH0uVrqYpqSFi6quU9UuVY0D1wAHu7tWAeN8h44F1rjbxwZst1gsFjtDv4iUtHBxfSgeZwBeJNkdwFkiUiUiE3Ec98+p6lqgVUQOcaPEPgbcXtRKW8qCG55ewfrtbX1dDUsfYRUX81T0dQU8RORm4ChgmIisAi4CjhKRmTimrRXAZwFUdZGI3Aq8CnQC56tql3uqz+NEntUAd7t/FkuCtzfv4sLbF3HbC6u5/fx393V1LEUkkbiyj+sxECgZ4aKqZwdsvjbL8RcDFwdsnwdML2DVLP2MLneJ2y079/RxTSzFxs7QLx4lbRazWEyQXNPD9jADDZtbrHhY4WIpOV56eyu/e2CJsfMn1/QwVoSlVLGhyEXDChdLyXHa5U9yyQNvGDu/NY3kz9f+toCTfvtYX1ej19hgseJRMj4Xi8VSuvzzxdV9XYWCYscV5rGai2XAkdRcbBcz0IhYtbVoWOFiKVlMdf427frAxZMtcXvzjWOFi6VkMTW4tGnXBy7W51I8rHCxlCxxY5qL89+Gow5c7L03jxUulpLF1OtvQ5EHLtblUjyscLGULOY1F/O8uHILS9e3FqEkSxisv614ZBUuIrKviNwtIneKyGQR+ZOIbBWR50RkarEqaRmY9AefyxlXPMVxvyn/+SH9Dau5mCeX5nI1cAVwI/AQcA8wGPgx8HuzVbMMdMx3ALaHGWhYf1vxyCVcGlT1P6p6M9Chqreow39whIzFYgxTZjHvrHb0OvAQGy9WNHIJl6jv82/S9lUWuC4WSwqm+34rWwYw9uYbJ5dwuVxE6gFU9Qpvo4jsBTxgsmIWizHNRb3/tocZaESKGMwx0MkqXFT1KqADQESqfNuXqupXzFbNUmqcccWTfPTaZ4tWnum+33YwA4/EDH07Rd84YUKR/yAi1TiOfcsA5sWVW3l8ycailWdKs/CcuVZxGXhYn0vxyBWKfCQwD3gcmC8iRxSlVhYL5vM/WbPYwCMxsOjjegwEwk6itOLeUnSMaS6ez8XI2S3lgB1XmCeXz+VR4CDgcGCOqtrZYJaiYdwsbjuYAUdyYGFvvmnCaC6fU9XdwBdMVkRErhOR9SKy0LftlyLyuoi8LCL/EpEmd/sEEdktIgvcvz/4fjNbRF4RkaUicql4+R4sZYc5n0vqf8vAwc5xKh5hhMuVrkP/csN1+RNwUtq2+4HpqjoDeAO4wLfvTVWd6f59zrf9SuA8YIr7l35OS5lgXHHpJz3Mns44W3bu6etqlAX95JaXBbkc+kdQJIe+a3LbnLbtPlXtdL8+A4zNdg4RGQU0qurT6vQcNwCnG6iupQiYm+fSv5y6n7txPrN+fH9fV6MsSEYK9pe7X7rk0lwk7X9f8ingbt/3iSLyoog8KiKHu9vGAKt8x6xyt1nKEOPzXPpJ//LQ6+v7ugplRz+59SVNWTj0ReT/AZ3ATe6mtUCLqs4Cvgb8VUQaCRaCgc+RiJwnIvNEZN6GDRtMVNvSS4zP0LddTN6U+4g/mZ2hb+sxECgZh34mROTjwHuBc1xTF6rarqqb3M/zgTeBvXE0Fb/pbCywJui8qnq1qs5R1TnNzc0mm2DpIVZzKT36yzWzAwvzVOQ6QFXbRGQi8EURmeD/jaqearBuiMhJwLeBI1V1l297M7BZVbtEZBKO436Zqm4WkVYROQR4FvgYcJnJOlrM0V86sv5EXJVISVjJe0a5a17lRE7h4vJv4FrgP0DcREVE5GbgKGCYiKwCLsKJDqsC7ncjip9xI8OOAH4kIp1AF4525QUDfB4n8qwGx0fj99NYyghTZjGP/tbNqCqmI+/LPSWXNYsVj7DCpU1VLzVZEVU9O2DztRmOvQ24LcO+ecD0AlbN0keYev8THUs/62BUk4kZTWFa4BeL/tGK0iascPmdiFwE3Ae0extV9QUjtbJYKIbm0r+6mGK0ptxli51EWTzCCpf9gY8Cx5A0i6n73dLHLFy9je1tHRw6eVhfV6Wg2KzI+RFXJWrYH1Lumkuy+uXdjnIgrHA5A5ikqnYacAny3sueAGDFz97TxzUpDCJOJ2DXc8mPYnT8ZS9c3Lt+83Nv86GDWpg5rqlvK9SPCZsV+SWgyWA9LJYE3tjblPO4v65EWYzm9BeHPsCfnlzedxUpALv3dNHe2dXX1chIWM1lBPC6iDxPqs/FaCiyZWAirupio8XyoxhaRX8SyNFI2LF1aTL1wnsYPaiapy44tq+rEkhY4XKR0VpYLD48zcVUP9ZfnbrF0CrKXnPxfY5Fy3e+jseabW19XYWMhBIubhoYi6UoJNY572+9v2GKoVWU/T3x1b86Fu3DivR/cmVF/m+uE4Q5xmLJB9PrnPcn046f4mgu5X3tFGfwMnpQNa1tnTmPt/ScXJrLYSJyR5b9AkwrYH0sFqu59JBiCM1yvyWqTqfVUB2jta2jr6vTr8klXE4LcQ4bnmwpKEnhYub8Zd4/ZsRqLuEQEQbVxti6ywoXk2QVLtbXYulL+qv5yhTFmedivAijePNcmhuqeG3N9j6uTf+mvGPxLP0Sz+diep5Lf6Mo81zKXLp4ZrHm+io2tLbnPN7Sc6xwsZQcyeSL5d2RFRvrc8mN59Bvbqiitb2T3XtKdxJiuZO3cBGRwSIyw0RlLOWB6U4sImY1l/4qtKzPJTeO5iI0N1QBsHGH1V5MEUq4iMgjItIoIkNwUsFcLyK/MVs1S6liun9JpH8pcxNMsbG5xULiai4A61tLdxJiuRNWcxmkqtuB9wPXq+ps4Dhz1bKUMsXqYIyv59LPsA793HgO/eGucLF+F3OEFS4VIjIKOBOwkyYHOMY7GDvPpUcU43KVfQSf59BPaC5WuJgirHD5EXAvsFRVn3fXrV9irlqWUsZ0p1+s3GL9DZsVOTeeQ39oXRUVEWHNVmsWM0XYxJX/UdW/e19UdRnwATNVspQ6xn0urkO/3AfJxcb6XHKjqghCNCKMHVzD25t39XWV+i1hhctCEVkHPA48BjypqtvMVctSyhjXXAybxcq8f8yIFS7h8J6vlqF1rLTCxRihzGKquhdwNvAK8F7gJRFZYLBelhKmaGYxo6X0P4phsip32eKvf3N9FZt3Fj571e0LVrNkXWvBz1tuhA1FHgu8GzgcmAUsAv5WyIqIyHUisl5EFvq2DRGR+0Vkift/sG/fBSKyVEQWi8iJvu2zReQVd9+lIlL+izaUGKY7seQ8F0OaSz8VW3YSZW6U5OBlUE2MrbsKL1y+fMsCjr/ksYKft9wI69BfCXwFuFtV36Wq71HVnxa4Ln8CTkrb9h3gQVWdAjzofkdEpgFnAfu5v7lCRLzFGa4EzgOmuH/p57T0EtPzT7zhQNlHJhUZO4kyN6pJn15TbYyde7ro6Ir3ca36J2GFyyzgBuDDIvK0iNwgIp8uZEVU9TFgc9rm04A/u5//DJzu236Lqrar6nJgKXCwGy7dqKpPq9Mz3eD7jaVAFG2ei80tlhfF0MjKXrigKZoLwLbdNjuyCcL6XF7C6dyvBx4CjgS+b7BeHiNUda1bh7XAcHf7GOBt33Gr3G1j3M/p27shIueJyDwRmbdhw4aCV7w/Y36EbDr9S/8kXoQBeL+4J650aap1hItNvW+GsD6XecDTwBnA68ARqjrBYL1yVilgm2bZ3n2j6tWqOkdV5zQ3Nxe0cn1FscxIpssxbRYr88F3RoqhVZgM3d22u4N/zF+V+8Be4L9ETbWVbrl2SSoThA1FPllV+2J4v05ERqnqWtfktd7dvgoY5ztuLLDG3T42YPuAwLEnmy/HvEO/OOX0F0Sce18MofnVWxdw+qxAY0Cv+e4/X+HOV9ay78gGpo8ZZKQMSI5Ah7jCZfNOq7mYIKzPZY+I/MYzI4nIr0XE3N1Pcgfwcffzx4HbfdvPEpEqEZmI47h/zjWdtYrIIW6U2Md8v+n3FMsebj4U2ZtEaT5abIuBUNREOb76m9T2Eok+i3D/xw6uMXZuLyx4u0EfiKomHPqD6xyzmMlnYCATVrhcB7Ti5BY7E9iO438pGCJyM47pbR8RWeUGDPwMOF5ElgDHu99R1UXArcCrwD3A+arqLczweeCPOE7+N4G7C1nPUqZYI/1yT1zp5+XV5uYC+y9TcVKzmCukocoxchy5tzkTcmWF0x21d5pzHnnpXwCG1Lmai4FwZEt4s9hkVfWne/lhoSdRqurZGXYdm+H4i4GLA7bPA6YXsGplQ7lHcXkUc4Z+l0EvuL/2cVUigS5BM+WZOndrW6exMmJR5/qYDg327kJNLEpVRaQsNZdyCNMPq7nsFpHDvC8i8m5gt5kqWXpK/zGLeeUYLcYpw2A/lmIWM1eMLxeb+Qtm0mRVEXG6o06DN95/iUSEIXWVRmbpm6YMZEtozeXzwJ9dP4vgzEf5RPpBIvL+EOdqU9W7QtfQEppimcW6jE+iLF5nWYzRPjjXLBbNeGivKIYw9u6FSc2logiai5L0uYCTDeLv81fx8w/MIBIpQjRMgSiH+UahhIuqLgAOEJFG9/v2DIdeg+NAz3aXjgCscDFA8TQXs+dPhiKbLQfMXjP/qTu64lQbki4JM6LJEb/736RwqYw6mktHl9l74u+cVm91DDDPr9jM3ElDjZVbaMohkjKrcBGRr2XYDoCqpi91fLeqfirHOW/Mp4KW8GiRslgUbT6NIb2iWI52f/13tnfRUB0zVximNRfn//Y2g2YxV3PpNKq5BIfrR8tIa4HyyI+Xy+fSkOMvBVX9SK4Cwxxj6Rn9RXNJJK4sgrDsNOnQ912nnXvMjfiT5Zm//2bNYq7mUqylToHvnrIvAI8v2WhUcBaaMrCK5TSL1arqt0Xkg/7FwnqCiByvqvf35hyW7PQbh77paDHfqK9YSQt3tpvrlJ15QWbHst7Zd7R30hVXIyP9WKQImkvaRXr3XsMA+N2DS3hh5Rb+8um5xsouJOUgXHJpLqeISAy4oABlXVuAc1iy0F/muRRzPZc9JudU+BqwfbdBzcWwMAanLZ5PZIch7SWW8LmYFPiaYharr0qOrxcanPNUaMrBLJZLc7kH2AjUiYjfiS+Aqmqj/2ARuSPDeQQoH29ZmVIsX8jfnn+bH51mMD2H4Wgx/2n3mHQe+zqATTvbjZVTlGgxoLGmgo079rC9rYNBtYX3H1X0gUO/zidcysnvUvaai6p+U1UHAXeqaqPvryFdsLgcDlwF/Drgb0eB625Jo1iayw1Pv1WUdCbFsFgVS3MxOZfCG4l/5ZYXjZWBQqObot6UbyKWcOgbFi4+GVJXmRQu5bSuYBnIltCTKF9P3yAiPw847hlgl6o+mvb3CLC4F/W0hKCYse8mU3R4o+In39xo5Pz+q2R2TkWSTTsMChdXHG8xnDq+0Y12M2XiK45ZLHm9AKpjyS4wWk7CpQxUl7DC5fiAbSenb1DVk1X14aATqOoR+VTMkj/FFC47DDqoq9wcU3e+vNZYGR4dRjWX4pjFioGiDHVzcW3YYaYtiUmURlPypL4jfm3lne1txldZLRT+WhbjPekJWYWLiHxeRF4B9hWRl31/y4GXwxQgIu8tREUtuSlWzi+AXe1dmQ8sICayAfg7/T39QXMpwoBbFVqG1gKwctNOo2UV0yyWzm0vBK8nc9mDS3h51VYzleoB/nf9/L++0HcVyUIuzeWvwPtwZt2/z/c3O4/5Kj/qefUs+WBac6nxzTA3OW+jmPNDjGbg9bVjk0GfS0N12CxOPUeB2soowxuqeGuTmQXDvOtl2lSZTRZn0sp+ff8bnPr7J43UqUeUgYKVy6G/TVVXAN/GaY73Vy8iLSHLKB9DZpljWqP3d5a7TAoX32cT80P85zeatNAvXAyZkgBahtQaO7eHqiII44fWGhMuHkaFi5aX4z4T5RCKHNbncifwX/f/g8Aywq+T8tke1MvSA4rhc2luqAKcdCbFwKQ5CWCjwU4/JRTZYDvSc5iZQgRGNFazvrXNyPk9f4fJUOQgTpg2IvG5DPzkQHnUM5RwUdX9VXWG+38KcDDwRNCxIvIT9/+P3N8+V6jKWrJjOoJEUfYd6WT9+dh1Bm+rrxmnX154U4R3maIRMStc3HKG1VfS2t5Je6cZgWxa0/OXMay+io2GBKVXhumsyOn8+swDEp9/ee9i/vrsSmPlF4oykC2hNZcUVPUF4KAMu58XkcuBeT2ulaVHFCPQxT+jua3DvPZiZm0P55zD6ivZ2GpQo3D/D6t3tD1TJjj/oGLnHkMCzJ182NxQxY72TiP3vhg+FwIc+ukJRb/7r1fMlV8gyiEUOZQnMC07cgQ4ENgQcNxFwBDgbKBTRGaqqnXoFwnTZjFVqPVNOrvxmbc49/BJhS8HZdSgatZuM2N+8WhuqGLxO60p66oXEq8DcIRLK5t27GHUIHNr0IPZHGaI0OwKyg2t7YwrsK/H0yr2dJpdOiDfW52y6JuhZyVf0q9QW0eXsSUdekpYzcWfCbkKx/dyWvpBqvpD9+O73O9WsBQR04t4gZMCxMNUBJQqTBxWxxeOmgwUfhZ90lxVRUeXss3Q6opJzcXs/BD/Xd9qcCKlAMManLYsWd9a8PPHi6C5eIEJ6Xz52CkZf+N/re54aY2JauVN+jhyy67SW00zrM/lh67g+A3wO1W9SVUzDSuvU9XFwHWFqqQlHKY1ZQUqK5KPTMxgLiYRf/CAmdG4Nwo35Xfx7seoJkdbeceQJua/779/eKmB8ycL8ExIX7v1pYKX4zXEZEoeCNZcvnr83onBTDr+9ucS3n211pGpZKK9IZRwEZHpIvIisBBYJCLzRWR6hsM/6P7/QCEqKCL7iMgC3992EfmKiPxARFb7tp/i+80FIrJURBaLyImFqEc5UIxoMUF44GtOsgVTy8J6rfCSChY6G4B3fk94bTDkd/E6gBENVUQE1rirHha+HDh0spMXduzgwpvdvMdKBGa3DAbgXQZWbfTuS7EmtqZz8MQhic9+K4D/NyWT29Kt1PQxTorHVpPm0B4S1ix2NfA1VR2vquOBr7vbgiioQ19VF6vqTFWdCcwGdgH/cndf4u1T1bsARGQacBawH3AScIWIlJYx0hDGrWLu+fca3kBVRYTdBh36gtDgChdTEymHGdZcvOtVEY0worGaNVsN+ZBUiUUjDG+oMjJy9s4oCJGIMHVUoxHTVbwImkt6VmQ/B44fnPjsj+zzX9JcA6pi+dm9d32fEa5wKVfNBajz5wxzE1HWpR/kOvSPwXHoHysiFxaikj6OBd5U1beyHHMacIuqtqvqcmApTuh0v6comov7btVURtltLDLJaUdCcynwi+NdpqTmYtYXIgKjm2pYu82c5iLiZC024T/y7od374fUxYxEvnn3xahwIfMkysbqGN97z1QgNUrRb4KKlIAzH5J1Guwmed1arj4XYJmIfF9EJrh/3wOWpx9UBIf+WcDNvu//6+Y6u05EvGHHGOBt3zGr3G39nmLMc/GoiRkULjgdWX21GbOYR1NtjAqDc10S5iSEUYOqjZnFnDKgsbrC6KJkXrc6uLbSSOBAUcxiqllThlS4mkmXbyJniuaSQ7YUK0DYq9NQV/teaThrQk8IK1w+BTQD/3T/hgGfzHDs9cAM4DYAEfmeiPxTRA7sTUVFpBI4FfCWW74SmAzMBNbirBkDwVpvt3suIueJyDwRmbdhQ7eo6rJk9x7zC6B4F7emMmrULAbJOTUF97m4b2ZEhKH1leaEC8kRv6O5tJkxW7mnbKyJGVlrJb3Gg2srjUQKeu1oL8L8qUx4C4Z1aSbhksssViyHvkNTbcz4ZOCeEmqei6puAb4U8tgFInKDqs4QkcOAE4Ff4QiD3ixQfTLwgqquc8tZ5+0QkWtw0tOAo6mM8/1uLNAtflBVr8b1G82ZM6f0ZySFwGS+L0h9yYxqLm45nlnMVLSY4PhdjJnFEpoLDK2rpL0zzq49XSmrHxakHJy5Fw1VFazYWPiMxX6HPjiJMo3kfHMLMrmcg0LWbIfRiDPe7spgFiuV1Sq9axUVYXhDFbsMvYu9oUcz9EPgtfQ9wJWqejtQ2ctzno3PJCYio3z7zsCJZAO4AzhLRKpEZCIwBejXKWgSo60izHPxOphYNGLUfCEiPs2lsC+O/yoVI52JSHIJXxMZBzwndWNNBdsNOHaTGphz86tjUTrjSmeB738xzGJkcehD0iyW4nPx3bJcwqXYZjHEmR5g9Jr1EFO5uleLyFXAccDPRaSKXggyEanFWbDMnwTzFyIyE+d+rvD2qeoiEbkVeBXoBM5X1dIT6wUkKkIXaihVShL/2WNRMbbuRiIUudIJ8jMWwy+OWWG5gdE++BzhSKLTMrEYlbdGSX1VjFYTZrG0Ksd8a91XFDAO07teRtdzIfsM+0jAffLXptSyrgiOQDR5zXqKKeFyJk4Y8K9UdaurZXyzpydT1V3A0LRtH81y/MXAxT0tr9yIRgS6iqS5uOO+aETMlec6XSuiEapjkYKb+/wdREUkYqwdiQUVJdlpmRsACJVRoaNLjaUoSWqtXlviQOGki3dpOuPm2gD5ay7+KMxc969YwidpqhRi0YjxpaF7QlbhIiKXkUXTU9VAP4wrDP7p+74Wx+luMUDUeMfl4HdWxqIRdnQajEzymd9MpWD3NIpOQ8vqfuAPT7nl+KKQTGgu7n/PXxBXiBp0DZgyw/o75q64JpY9NlVGEMm2JZ+J1Hplf1aKtc5KwlQJJStccpmq5gHzgWqcZJVL3L+ZJP0qlj7G64hzPfiFLMukKu4/q6MhFdq273PQRgVT76UXKCAivgFA4QtzRvnJNegLXYY/pBqSgrLQQt9/X0wNlHItc5wULikVM16vfPEHWVREpWTq5Ser5qKqfwYQkU8AR6tqh/v9D8B9xmtnCUWQKm+C1E4/YrYDcD87moUhzUWc85sWyqY1F6+MaMJfUNhz+0OqIRmcYFJzMfZsEZy40iNoEOAXernaXDSzmPtfBGKR0tRcwvpcRuNkRN7sfq93t1lKgGJFi/k7fcehbzZaDJy2FTzzQFr0j2mhLGL2Hnmj8YqUjrHwGY+8e29KC/ObXU09W7k0l6BBQIrQKxHHuT9Y5LkVm3Mc3TeEFS4/A14UkYfd70cCPzBSI0veRA2ZKQJx38yKqEHNxW+2EnPmN0+jMC2UTQsXcDqZYvhCwCfECm4WS2Li2Xpl1Tbue3UdY5oyJ/eMBFxD/+CmGCmWwuDXXEqVsJMorxeRu0lOgvyOqr5jrlqWfIiKuTDXTJh0hINvlBwtfOdfLPOehz8U2cQAwBPGpsyj6R2ZqQAS0xrCJQ+8AcDqLGl4AjUX3/5SixZL3VYaC5l5hE25LzhzVg7wJkSKyIBIBlkOmA9zTZLiCzHl0E8PFTb0xjqOdvPmREdz8SK5zJnFgmaXF+b8SRMMJOe5FLqc1JDfwg9cwtQ3SHCmR7GVBp4fLClMSqduDmEnNl6Bk4zybPd7K3C5kRpZ8iY52jKb8C+lTIMhwn67eETMjpCj7jwX0zmhTAZdKOk+l+JoLoV2IqdoCAaerVDCJcAKkBLFlqNeRQtFdouJCJy43wigdCLZPMIKl7mqej7QBolcY71N52IpEJEi+lz8ochmo6y8sNdISobagpYg5qO4nHIkcP5EofCW7k2UYfg5MHXNUqPFCn+dwgjDZDi3rzJ5zHMpFl71BGG2uw5NuWouHe6CWwogIs1AaVxlS8JUZfLh6ubULUL6FzATzRWUiNCU6Q1SQ5FNXDN1CzE9z8XDVJ60lGgxA89ymHMGJ64Mf47ihSInw8O9Oper5nIpzuqPw0XkYuAJ4KfGamXJC++RKo7PJWl3N1me3wRjKkKnGPNPoFjRYgbL8EbJXqRgQlCaE2ImhHCY+npmsUzRYjnnufSwbvniz7htMm9dbwgbLXaTiMzHWQlSgNNV9TWjNbPkjVGfS9r3qMFoMf8I1ojmkuJzMR8MIUiwuaVQuKc053PxHPoOpoRYPjm8ekIYs3Euh36paAf+GfrFSv+UL6GEi4j8xU0U+XrANktf4z5TRdFcvJxfEcNJEt3/Jn07KT4Xw34KU5Fc4Dn0xWC0mPO/e+JKcyHiJiZRhhkMeYOATGax3DP0i+TQT4ZZpE2eLR3CmsX2839x/S+zC18dS0/wHjOzPpfu0WKQdCwWmkS0mIFJjik+HYPrrCTLU7PRYl4WacMj2KTmYj79i4nglNNn5V7t3FtpMnUlyvBmsWIRqLmUSPYAj6zCRUQuEJFWYIaIbBeRVvf7euD2otTQEppiznMxFY4K6fNcTM6gl0Abe6GJa3C23ULhhSKbKiP9ylQYuvemO/GGECuABoX152MW62mtN+1o75G2JgRrW6VAVuGiqj9V1Qbgl6raqKoN7t9QVb2gSHW05MB7KU2adtLPbMo04pSVTC5oxueSPF8xTAqqZjUXMBuRlphEmUj9YygU2fe5w5AQzoVfC1iyrpUJ37mTJetbE/tzptzvwSWJx5XZP3mAr//9pdC/8a/nUlHO0WKqeoGIDBaRg0XkCO/PdOUs+VFMn4v3QJucgwJutJjBrMjFSvppOnGlyTIS1n3fHCco/POWEpVl4LkK8xwlMkurcs9CJ8PVP19YndhvJNDAFVi3L1gT+jf+IItS9bmEdeifC3wZGAssAA4BngaOMVYzS2iSocgmR9+p373Rq5ERZppZzKTjuBgmhbhqcnRpKLeYiOGINLr7XEzOpzGy7k2IY/yCMxrwbOR8Tny7wwa79KSpZe9z8fFl4CDgLVU9GpgFbDBWK0teeA9acTSX5Mx5MDcp0K+5lOXcEB+qBHZUhTy/0xbzjnYwZ34zHYoc5pT+rMiRAH9cPvUKex96MoHXr02Wpc/FR5uqtgGISJWqvg7sY65alp5g1ueSHi1mzqEPGPW5+JvidZRnX/NsYcvwEVfDucUUZ4a+sWzFxfO5uE0wM2gJ0Yl717C1rTNRl5SU+zkd+sn9e0K+Gz25jv5kohWGNMneEla4rBKRJuDfwP0icjsQ3kBoMYr3QBfToVfphvCaiRZLtiMaiRj0uSTnhmzc0W5sjoKqGo0Wg/T1XMxEi3VLXGlA6MeMPle5j/E0l1/euzihuXT2UHPp6Ax37M72ztDn9EicWcxpkr0l7Az9M9yPP3AXDBsE3GOsVmmIyAqcTMxdQKeqzhGRIcDfgAnACuBMN6EmInIB8Gn3+C+p6r3FqmtfYjYrcur3qgqnE2jvNFSmz3lsMreY92KCM9Ksqij8Co5qWHPxMKe5OP89bTIZzFHYex9XpTIaob0zbmiyae5zeuHKB7Y0JTS1fHwu/vckrOZy6M8eCnVcUDnFMu32hFzzXIak/wGv4OQWqy9KDZMcraozVXWO+/07wIOqOgV40P2OiEwDzsKZ+HkScIU76bPfUlyfi/O/KmZOuPhbERGzPpeIT7gYW0IANdoBrN66m9ff2W4wWswziznfTQUOqELMHbQUXCsinM9FRJg4rI4xg2sDzWL5zHMJK1x6RtJUmQyuKS3hkktzmY/rXw3Yp8CkgtcoPKcBR7mf/ww8Anzb3X6LqrYDy0VkKXAwTnRbv8R79osxcvFGr94If48JzUXTFiUzGJXk11xMrttuOv/TojXbjQVZ+EfJYDYU2Zs/VWitCJLt+N1ZM7MeV1URYU9nV6BDPx/rgKnnCVLXc0lokiXmc8kqXFR1YrEqkgMF7hMRBa5S1auBEaq6FkBV14rIcPfYMcAzvt+ucrelICLnAecBtLS0mKx70TCawiSjWazLSHmeSSJWIcY0Cn8YJ8Did1qZO2lowctJnaFvpi0fPWS8sYi0TD4XE+XEDKbj8TSQk6aPzHpcZYVjmvMejZRosVyLhfleFFPBLpC6nkuphiKHnecSOGFSVR8rbHUy8m5VXeMKkPtF5PUsx2bSslI3OALqaoA5c+aU1l3pIUXRXDyzmKu5tHeYNYtVRqMF144yaS5f//tLPPHtwk/dcnKLmQvdrqyIUFdVYT5azFtuwZiGpL5AEXPRYhLYRSSpqojQ3hFPDHA68/C5+DG5eF8ygq90Q5FDCRfgm77P1ThmpvkUaRKlqq5x/68XkX+55a8TkVGu1jIKJ98ZOJrKON/PxzJAItuMquFp8rnSoEPfS8TolWPE9EbqqM8kccU3CjZ0vUzmFkuoLs6/SEQQMTOJ0tNczFwn53+uW15VEWV3R7BZLB+fi0nNxX9LvEHF1t0dxsrrCWHTv7zP93c8MB1YZ7ZqDiJSJyIN3mfgBGAhcAfwcfewj5NMpHkHcJaIVInIRGAK8Fwx6tpXeKOYYiauNG8Wc/5XVkTY0xUvaJiw/0yeRgE9ywsVrkB3Bn1EjDhdVT3bu9losYhvtrmJKL64KrEKb/6UOYd+rlnzjlmsK2FGy7RwWBCpmZ2LEL3p87lc8M9XeH7FZmNl5ktYzSWdVTgCphiMAP7lPhAVwF9V9R4ReR64VUQ+DawEPgigqotE5FbgVaATOF9VzfSAJYL3nBVzmeNiRYt5QsxEmLBIaocZCTXUyh+vPbFoxIh2GVdnNrnxaDHftopIxNAkSjGWlSGoHUFUudqyJxx6ukJmGAH5wKs9G6Mn25Kqfb/w1hYOmjCkR+csNGF9LpeRfEciwEzgJUN1SkFVlwEHBGzfhLMyZtBvLgYuNly1kqOoocgmo8VIdgBe9NCezsIJl5SsyNHki3nSftkdvT3FE2BVsYgRYRzX1Oy4hR71py8WBq7mYqAccYWLiZx13usRyWEX8xz66c92ZTS3QPWbj8NoLufeMC/nMRkKApx74pkSwdz6Sj0hrObivwKdwM2q+qSB+lh6QGKeSxFsvB4mzWL+kaJJBy+kvpiRHOaSnpIUyIX3HyWd1Ek/WMHLcP/7L080WvgQ8bjra4sZEFxAaLun59BPaC7u9liebS6Wz6W+OtmN5zLbFZOwM/T/bLoilt5jMjrFw4u0STj0jUSLJbPJVhrQkPxXqd63gJSpbAMJzaUimrGM+W9tZun6HXzooPzC4v3+kGhEiEWFtgIL/KAoq4pIxIimHBFzyUqV3M58cO7Tnq6k5uK1P1YRwhSY4nMxb6YWEWpjUd/2MhMuIvJe4MfAePc3AqiqNhqsmyUknipuyrkO3R/aiogQEXMdsj9aDMyY30SgwTfqMzWjeuKwOsAdEWe4Rx+40pnje/iUZkY31YQ+tzdS9TrN6ooobR0FFi7u/3SzWKETpcbdwIdYNGJk1O+dPxeVFRHaO7pod+vgPXthBGpvosWmDA+f9MSfNcFv5isls1hYF+ZvcSKyhvpWpLSCpUTw+v02A1pEOt67KSLuSNywWSzh0C9cOf7ze+Y9KLwAmz6mkZYhtRziTsysikVyanrf/Ed+rsx0P0JVLFrw5yBoMGzCL5KIeosa0lw0rObiRCh6iSe9AVRlnvXKV7jkk3o/PWuCRymZxcIKl7eBhVpKOpclgXdTjGouAdtMOai9FPKQ9LmYKEeQlJGsicma/tFoNrOYR76mlPTOpDrmjLoLizdKTl6rWLTwfhHH5+IEJpgKRc41gRKcAU1HlybMi55GG8YslhqKnF8b8rmeQdoklJbmEtah/y3gLhF5FGj3Nqrqb4zUypIX3gPd0aV0xdXIxMCgYYUJB7VH0reTjBYrHMFvoAnh4n/5s5nFPHo6Yvf8OtWxqAGfi1dGcpuJya2JdWkMBAuA58vLfZwXleilwt+T0FwibOsKP1Exb80lj3ufHOdLhu19T1jN5WJgF87s/Abfn6UkSD5QxrSXgAloYUbivaUymt2hH49rj6Pk0juaQvtc0mdVVFXk1vTy7bC7+VxikYKbxfx5rDyMCDGSDn1Ti6qFEy5Ot9ja5giXbe7M99qqipz+rJRQ5DzvZT7CKLPmUjrCJazmMkRVTzBakwHE5Q8vZe223fzk9P0Lcj7/89TeEae2siCnTS0jYAJaZYiReI/KUk2ZoQ+ZTQzvv/IpXl61lWU/fU8e50/9/tz/O5azr37GSJhwquYSzelzyV+4OP8TmosRh37SeexhpBxVRCLEImYmm6pqqHBz75nb0Za6iFddpePPisc141yZ3pjF8rqeGX0ueRVplLCaywMiYoVLgfjlvYu58ZmVBTufkkz9UejRZKKMgIl03nwAE3SLFsvg0F/w9tYev1BeW4Y3VNPcUGXExOfvgxwfVfb7k2+qG2+k6rWlptJEp+/893dk1ZVmAgciEXOhyI7PJTcJzSVthcjaSmcsHlZbz1cT3rWnK/S91wA/GJSW5hJWuJwP3CMiu0Vku4i0ish2kxWz5Ed1zFyWYkidtOURxsxz1ytr2bJzT4/KAsfMA7B7j5l5Lh6VFdFE6GnByklzIIe5Xss37mTiBXfx3PJwOaLUPZ345tLsNhQtlqq5RAouxDyHfmWI69QTvAwAuUhoLu2p/pX6Kucd252l3f5na/ee/K5PZ1xDtztd4P/hI7OB8EsrF4OwiSsbVDWiqjU2FNkMj76xgasefbNHv1XVRCcc5uFUVVZu2pV3GRDkc8n8Aq3dtpsv3PQC/3vzC3mV5ZTj/G+sjgHQ2lb4jK/+jr8yasBBTYBZLGQZdy9cG7oMSGpINZXRwGixn971Guf++flQ58xUBuk+FwPzaUSgJhbNu2MOd/78HPqtaWYxbwCXTbj4CXucn10h252eTPSk6SMZXBsraMh+b8nqcxGRfVX1dRE5MGi/qubfa1gC+fh1TuLmzx45Oe/fKskXIswLf9OzK/nevxdy9D7NXP/Jg0OXAWlmsViEXTs7A4+HpP9g5eZ8BVnyszfJcXsBhUvmyLdCj8RTr1dlRYTNO/cERvQ11cbYuivZxnhIs1B3n0sksFO76rFledY+SaDmYihwQESorYywdlvhBxMa1iwWC/a5eOaybILPb9YKKyggeT13tncypC630zTdHApml6foCbkc+l/DWa3x1wH7lCKt51LO/PaBN7hjwRoe+sZRoY7XkLOIU3/jjFghnOay4O2tADy8eENe5UCaQz+a3XzhOTTzNQMrmtAq6ioriEj3UWQ62ZysmUh/Md/csJMtO/cwOMTLHQZVTdGOHn3Dud5/emoFnz5sYsqx6XMcwk6o6x4tVniNwsN/dWsMRIvh5harraro0ag/9+nDPSNVGVbD9Mxl2a6v/7bt3pP9mfVT405+Da25BGwrNeGS1Symque5/48O+LOCJQS/fWAJyzbuDH18T2zNfrNYmI7Fn6wxfBnB53n9ndaMIzmvLj3xMXodfyQi1FdVsD1tIaStu/akOH3zcZ6mL3wGsNn1C33hpsIp4wopPfL67W0ArNzU/XlID0MN25x4mlpRUxnN2kH1xFHuz2PlUR3CdHXD0yt4YsnG8OXgCMnaWJRdeXTMYQnt0I8Fvx/ee2PCLFbn5rjb0R6u3UHaZGU0YiyFUU/I2suIyEEiMtL3/WMicruIXCoipbFoQD/j3kXv5P0bxQkNBfjYdWHWRetBB0P3p/kxdyT+83uCV51OCpf8yks/vLEmxnaf5rKzvZOZP7qfyd+9K7GtJ4EM/o7GM1OtCOj4e0xaZ3bkPsOBpO3eT/ooOaxZzLst3oC8ttLx62QSIj0JHU/364CTZqa9s3tk2xNLNvLPF1YBcOHti/jItc+GLsfL/ZVLQPYUJVwocqalHSpDmMX85NMGL4FqeKHqTQ1I9YGWjeYCXAXsARCRI4CfATcA23DXn7f0nKBO98u3LOjBiZJmsTDs6UlESUA4qneWTRmiwbwHPZ+cSd55/X1AfVVFilks6KVtz8ORGVQdr+MsZCSn045kQy587zQAmhuq0uqj3YRBeLOY89/rNGsrszud3968O9R5g8rw35OaWLAZ9iPXPsvXbn2JCd+5M+9yvNxftZVmHPrpPrBMVPryzVX6tPx37zUMCK+R5GOe9DSXl1yTdS4CNRdDUXY9JZdwiaqqFxP5IeBqVb1NVb8P7GW2av2fQs1CVhzfRFjGD63tURnQfZ4LZDY1XPvEciC8iSeV5FkbqisSqTjS6+DRI83Fvz6JgSWCnUl7ye+Da53It3ThGDTZLux8hXSfizcXY1cG88qJv32Mt/MOsOg+SvaE2M6QZpwwxF07Ym1llM64GkrHE0ZzSXaLQ+sd/1tEYFi9MyjI7tBPfs5Hc/GE9a/ueyPU8UHvY1n5XICoiHi91rHAQ759PV0ieUASZOYoVOI/VU1JHZ8Lr2MY0ViV40h/Gc7/lA7GjfsP6gg7uuI8+Pr6jPvDlOVRV1XBTp+5IMjkk5/PpTveKo6F7CyVtKWBoxGqKiIpbQEC82iFvWTxtI6/zr0n/o4tXUPOxwcIvuvla0yTKyi37i5kFJ8jjGtcAVlo7UXdgIFc+DUXL3IrrkkBkE1o+P15+dR/7GBnmYVJzXWhjg96H71szqVCLuFyM/CoiNwO7AYeBxCRvXBMY5aQBI2IC/kgpKxGl2P07dUln/QUQSlAhjdUA8Ev20OuYIGezBpOnY9QV1WREhYadC17Nnk0WYgXRZTJ5PHw4vX86cnleZ09aKRcV1XBrvbCaS7p5pGamPMcZBPG+Y5ug2boN7k5hvKdIJurHBEnzQrQTQgX6vy58Ptc/GHB+YTFV0QkL4f+ns44cycOSWhHuQh6H03M1eoNuaLFLga+DvwJOMyXcj8CfNFs1foXmUb3QeTtAAeivqcsV4iol7cpr0R5AR3MZWfPApKjLj9+M51/tceeUF9ZkRJFE5R3Ki/NJeD6vmd/J24l08v9yeuf5wf/eTV0GZBcttdPbWW0m3bktWfckOR1DBvVlT6ZLkhzSfff5B+J1X0C7RBPuOxK7WjT/UlOHUO2xQ1B9/yHhXbqh3foJ7tF//NQWxklFpVEIsvAMtym5huUcNqsMTTVxti2K5wmGPQ+lptZDFV9RlX/pao7fdveKNYEShEZJyIPi8hrIrJIRL7sbv+BiKwWkQXu3ym+31wgIktFZLGInFiMeuYiaLSdySyWr1POG5F9zp2AubM9+0PtjZR7m4XVWzHxiaXdw02vfjw5aa8iz/kn6ZPd6tN8LkHXsifmLH9bTtxvJOOH1rL3iPCrAebCvy6NR0N1rFvOKq89+4xIJr0Ie28SPhf3TU74XHwdW7rVLVvnGESw5uKYxbbsStVcagIi4cKO4L3cYrUhzWIrN+1KhJCHIWwost8sNqgmlvgsIgyqiYW6fvVVzjObS7DWVUb59GETOXLvZppqKtm6O1x7gt7HmoCBS1+S/4SH4tMJfF1VpwKHAOeLyDR33yWqOtP9uwvA3XcWsB9wEnCFiIQPpTJE0Eg0UweydP2OvM7trTnvLUyVa2TqletleA1VRoBT12PZhlQb/p7OeCJM2alP/iPQdLPYzj1diboGCeX0Ti7/8oSWIbUFn7yXfr0aq7vP2fFGm8N9PrCws9/TfS6eP83v0E/XXMKOjj2COjJvoumqtOCAIOGSawKsh5dbLNGGHM/xEb98mHf99MFQ54bwDn3/PLB/L1idsq8xh3DxrtWgmhidcWVnjmc/rslgkkG14QQX+LXBVG0yXZPsS0peuKjqWk9LUtVW4DVgTJafnAbcoqrtqrocWAqEy3FikFxOaL854YWVW/I6tzfS90wi2TSXn979Glc8ksxhFjatSoa1iTh95mgAnnozqb34HdTvnTGK7W0deZn60o/0EgZ6NvggB3hPbP/p3UyYnFb5pIJPT7kP3efsQFJT3du3amXoyXTuf6+c+oDJeOnPXk81F79JyfOLXPrQ0pSQW397j5vqzOsJK1y8EPSkbyP37/LR8oPuRy7SheWgmuymK0/Ye76aXNe6y1enQTUx2jriec1F8rdnSH0luzu6jIRx94SSFy5+RGQCMAvwZmb9r4i8LCLXichgd9sYnGWZPVaRXRgVhSDh4ndSx+PK2Qe3EI0I67e3J2Zzh8GN4EyaEzoyv5RXPZqaY2rjjvw65fR3c0Sj49T/8DXJyXL+zKwtQ2rzSmsB3dOm1Fc5pglPaAZdy807w3eYmeRc+nyaIHblodkEWMVorI5l1FxGDqrmtR+dxHtnjGLTjnbC4Altr+P3zFX+Ts3T+Ia6Hd7mPLW8pHaUxK8BZBJkQ+ucAVPYpKOeZjE44c/JXE+/byHswGVPVzxl3koYrvrobK766Gzu/NJhADnNYl5VvDZszXGt43FN+Esba7rfu1zl+O+J5wfL9/6aomyEi4jUA7cBX1HV7cCVwGRgJrCWZP6zoLFJt6dPRM4TkXkiMm/DhvxzbOVLYIfoewg6uuJUVUQYVl/J7x9eysH/92DKDPSsKCnmhFw+l5Q6hBzxB6UAgeTkLz9+jWzckNq8yvHwF5MepeOPrqqICI3VFXmZxTKthTGkrjLneZasC2+ydCYFppnFaiq6aYtXutmwY9EINZVRmhuqQgv99EmU9VUVxKKSMrHVM4t9+bgpHDCuiQ2t4QSXvx1ARodFJn+Y57sIrbm4ocjeqN/fMb+yaluKEPGbq8KaMrft7qDB50MJw+Tmek7cbyT7jR4EhBEuTh0TQj6HmSquySSmnn8nffARWE5QkEVd4SP4ekNZCBcRieEIlptU9Z8AqrpOVbtUNQ5cQ9L0tQoY5/v5WGBN+jlV9WpVnaOqc5qbm802ADjzqqd5ZtkmHlmcDNH1PwQdXUosKinmsdARQ27obtKZG96pF3aE7BEU/ZSO32w10tVsNuZRTnqrvYlsm9wO139d/vLpuQypq8yYJSDw/AGjPnD8CLv2dGWdWb1qS/gJiEEp3hurY7S2dSbasHrrbv7zUurjOay+ih3tnaHMG+mTKL2Rv//Z8jSXaEQY3lCVv3AJSDXix9Nc3t68i+W+OTRH7t2csj8XnsO9tjJKZTSS0EYfX7KB9/3+CW58NrnAnn9ZgbADl+1tnSkO+mz87bxD+OS7J3QbPDXVxLJqI96TmRCQOQSRlwkaksJlawi/SaDm4paZz7tgkpIXLuJc+WuB11T1N77to3yHnQEsdD/fAZwlIlUiMhGYAoRJuGWUlZt3cdbVz6TM//C/FJ3xOBXRSGLuSD7k43NJJ+yDGJTsEeBU1+fiCRFINYslHvg8zG/p0WJeOOimnU6n6AnoW847hHdNHsrgusq8RmtBqTMgaTZKvyb3+NZWyaecoHQjicmHbgf17p8l5yUfPsXpjJvd9oYRyJ4cT7G911WmPVuucBFHuLyzvS2/cPcM1+sX/zMDSD5vx/7m0ZT9U0c70W+hzWJuYIqIMLgulrjW3pINr65JTq0b5FvLO6xwad3dQWPIycZzJw3lovft12370Poqtrd1ZhyAxBOai6d9ZfPPOP89s9jgRAReeOHi14w9E/VTAdGbfUHJCxfg3cBHgWPSwo5/ISKviMjLwNHAVwFUdRFwK/AqcA9wvqqWhocLuOHptxKfN7S2E48rqupqLpFEx9ITmkLYqj28LMphO/1MHfLwhmqO3Xd4ymSzDp/m4mkd+WguTjnJgoamCSgvIMEzJwyrr2J9a54+KrqPxL02bE67Juf/9cXE53x8VOo5w3yMGuR0AGu3pda3OhZJmJE87TVMm4LMI4NrU817npYUiQjTRjeydVcHKzbtYmd7J9/6x0s5/QKZrGLTRjnCwyvL7weZ1FyXGIl/+7ZXQgVC+Cc5Dq6tTJiNvQmyfo31Szcn70nYAdI729tCay6ZGOnev/Xbg59n7z0ZUpfbf9KV0Cjdc7vC4Z1tufO/Ba3nMsadGlAqmkvJp3BR1ScItvZmdEi4kz8vNlapArFs404+fv1ziRnBsYh0m4QWZp0SL8qmsbqC6liEdSGCAeLqqOEbdoTrlIPCUT0aqit4Y33yJfKHWI9orCYisGZr+ISJ6aPqptpKIgIvrtzCxw+dkNjudWZjB9fw1NKNodfC0YAXE/xmhdSOY7/Rjby8yhk15ycku5vFRg1yOoCVm3cxfcygxHZ/6PEYd1Lq6q1tzB6fo4SAEeyQukpefye5CrnXEUVFmD1pKADPLNtEe0cXt85bRV1VReAoPb2MIB8VdNcc/mf2WH7+gRkpC6K9tXkXk5uzzyFSTQp8Z8DQztL1rfzov87k1c4upa2ji7n/lxp+fN0Ty+nqUo6dOjzj/W/r6GLXnq5ui7TlS0IAbG+jJSBHn3eta2LOu5jteUnOUUq2ORYV1mwLM6hw8Dc3EhEOGDso1PtfDMpBc+m3rNy0i8eXbOSB19YBTu6pdOGyI4T/xIuuEhFGNFazLm1U1drWwYTv3MnVjyVDkPd0xhk7uIZVW8J1+tnmuTS4fgQPzyx2wLgmYtEIowbV8HbIcoKIRoS4wr8XrOGh19cltns+ibGDa9m5pyv8HIEM2z2zQrpWMbi2knFDapg6qjE/31HApL1RTU4Z2daN8UagYfw76T4XgMF1sZTRqxcAUREVJg2ro64yyhvrWqmKhUs+GZRqBDILl/bOeLdOPJd2BKmJPscNqWH1ll3cOm9Vsh1xZfXW3d3u8+NLNnLuDfN4PMvaMZ6T3Jv421M8zeWdDB14Utg7Za3NooW8uHIrkMy8HIkIIwdVszbEQCw9StBjWH3mYJDlG3dyxhVP5uU37A1WuPQRM8YOYnXaQxSLOjZxP2HMVp7mAjCiobrbyMU7x//dlVx35fpPHETLkNrQSxBn01yaap3wWk9j8cxiXz1uCuBoFvlk4s3mDfBPMPWcrV76mbCCMpMPYUxTDdWxSLdJrNt2dzBhaB0jG6tYvTU/81t6GcPq/BMlg621dVUVDKmrDNWeoHT4Y5pq2bqrI+FI9+YgxaIRRIQxg2tYvWV3qESMkDkAojoWpa4y2u0Z9TvbD53saEqZzEjpbfHaMXZwLRt37OFq3/LMXfF4SpqjdN7JMuL35sz4fYM9wRuArMtQll+THD2ohjUBz8sVjyzlqaUbOfuaZ4DUjACjMvwmnfQoQQ8n0jD4Wt82fxUvrtzKHS91i28yghUufcSpB4zuFkVTWdFdc1m+MXfoq3+EPLzRMSf4Sc9rNnfiEI7edzgtQ2pZtWV3qFn62fy/4wbXEldY674UHa65yhuRjRtSy+J1reGj2AIc4Q+7y0RHRHjfAU4QwbvcjispXMIKyuBQ5EhEmDC0jhVpWYO3t3XQWBNjcnM9r63dnkfeL+2m6flNnNnMF+NCCuSkiS953hY3/HvlJuf3F96+CEim4RnTVMPqrbsTQiVXVFpyYNG9Yx8xqLqbybPN53u51M0/l/5MBpeTzP01YWj37MAi0i2rxQ9PTZrzsuWXO84NNggKnc+HxuoKGqorMg7K/Fre6Kbu1wbgF/cs5sN/TM4L81/V0YOquw06gwjSWMHRXDbtaA98Rj1tMj1xqimscCkiY3wq+eTh3e3PAkxy7dLvmeEEw1320NJwJ3dfylHuy+4XGOmpRPYkEiXWsqczHurFD0pe6DHWTbj4ttu5e9FJFa5waRlSS2ubs3pkGOLaPcHghKG1VEYjbNyxB1VNSU3uzaVZvjGkcMkwEgeYOKyOpRt2MOE7d/Kb+xYDTnqbxupYwqTyh0ffDPhlQDl0f/kBrv34HCDV/HbXlw5POWbs4Nq8NBf/9fLW61m5OVVIemlNxgyuYcm6HXz3X68AuTWXIOexx8Shdd1W7/RrZENqK6mISGg/oFdGUOr59dvbUmbkz2ppSkn2GcYs6kVU9hQRYXJzPW9uCB70+e/H6KYa1re258yCscPX2Y8fWseabbtzLjSW1FhTb0rLUGegF7SiaizqCpcizeC3wqWIeLmjvnDU5IS5wM+uPV0Mqavkoa8fySVnzgQcu2y2JIZJX4jDxGH1tHfGWeOz9aZPMvvkuycCMGmY8wIvWd+as+7ZOuTESNkdzb3w1hbnWPfgA8Y1AeFTvXf5JpZ5iAhD6yvZ0NpOR1ecWCT56DZWxxjTVMNinxM7a1sS5+y+b1JzHW+5I/7fP7yU19Y65/z7vLcTYdfz3fblIh4PDjDwtFP/qHZMmi+guaGK5Rt3puRoC2xLwOx5T7i8kTbhs8LtXCYNq08Z5efUKLPc+wnDHOHiH8z4j4u4c2vS/ViBxfii6yYOq0vcn2P2Hc6H57aw+J3WlNQoXz1u70SABMAv712cs4zeZugGZ2Jlpvx/6ntRprsTL296ZmX3/T78wmefkQ2o5s4vqBk0l/3dIJGFq7el/ySRMDWM/6sQWOFSRCYNc7SSuqoKqiqizB4/OGW/l+RuUnM9lRWRxMz0bCPY9BDhye6Iz59M0m/2eM/+ozjVNSt5s45fCXgQu5VDajl+Rg2qoSIiCeFyz6J3gGQI8RxfO8OY4DItR9vcUMWGHe10dimxitQD9hnZwOvv5BaS4BeU3QvxUs149fCc1ROG1TGsvopj9h0e3vwWMEMfHDOiSGoOufRlqo/cx5nz8kd3Nc9MJMNZk+U0VMeY3FzHy6u2plxvb0G0/ccOSjlHriSZmcyI4FyXto4472xvS9Thy66vzWPfUY0sWhPiGfM59Ktj0YTAbayuYOrIBra3dfKP+cmZ+bWV0UDzWTp+E1FvzWIAk4fXZdRI/JqLtyzyz+9J+jqDMnofMSU5iXufkQ0AOZ9l776mP197DXf6jiDhsn23I1xW5RG52RuscCkiZx88jg/NGccH54wFknMePNrTNIzrPnEQAG8FqLjpeB2lZ1bzj3z8mstIX5mDamOMH1rLy2+HefFTy/ETjQgThtUlUqMc4c7MnuhqRnVVFZwztwUIF4Pvz7fkZ8LQOpaua2VPVzzRUXrs6wqXUJ1YhugngDPde5M41m33j05zbPv7jW7kzQ07Q62P3hnXhLbgZ3BdJeOH1HKjO6I9aMLgFKcuwNH7DGf80NpAs5qfTI7dA8Y1seDtbVznW+DMM4t481M8Vm7elUM7dv4HVWWq2xn+Y/4qKqMRPnP4RA6dPCzlmP3HDGLp+h0hotJS74kXulxTGWXKCKecm59LagH11RXdhHJQGf62FUJzmTrSuX6LVgdpyklN0l83b56PP6P38IYqPjh7bMJ3CDB+SC2VFZGcWnim+x6LRpg6qrHbgHFHe2diMvCra7aHzobeG6xwKSLNDVX8/H9mJGbhp0eufNZdj8VjL/flChqFeKQ/IsPqKxlWX8miNcmH0y9c0tcrObBlMPPe2pJzxna2Dhlg+ujk6LS1rYMhdZUpI92Tpzs+pIUhOv94gFkMnI59zbY23tq0K9FRJvc5o/H3XPpEzvNn6yyH1lexr9thQjITc2N1LFFOV1xTrm8muuKZF6dasSmp/Zx9cEvgMbPGNbE41wg2g3lkVstgNu5o5yd3vpbYNnawYy5LH73v7uhKzOMJIkg78pjuM8O0d3ZRHZByf8bYQcQ1+3MMnuaSLGOqKwQjIhw0YUjKseOG1LCv28n/5PTpiSW7F6/rfr38wqUQmsusliYgOHt5eqf/mzMPAEikxfFPMF7f2k5DdeqkzopohGmjGhNhyplI+MECevDpoxtZtDopQOJx5at/W5CY+b+jvZM3QpjCe4sVLkUkfXTq+SK+eMxerPjZe1JmuYMzwp02qjFr/H76hEAR4YCxTby8amvimDafWey9M0an/H7uxCFs3NHOmxuya0fZOmRwOpm129rYuKOd1rbOhEnPY447Os+VmiI935Ifz4y4cvOubp3YMfsOz3relDLc/5km3PlH9p/9y3wgmUNtzgSnDs8u35SznC7VjAulffPEfRKfg9ZAAcectHZbW9b0Jpk6/rkTUzvj9x84JkVr9Xx+fzvvEMCZVJlvGeCYr06YNoL7Xl1HXAkULge2DCYi8NSb2a9Z+mJeXif+5oYdRCPCJ3wTaG/41NzE548cMp67vnQ4IvD4G92fL3+i09oM1zofmmormdxcF+h7889zgaSAfNX13aWvRVRT2b0LPmTSUF5atTWrLyxo8qzH/mMG0dremTBT/99dr3H/q+tSjnlksflkvVa4GCRd9UxP933y9JF895R9U16adA6fMowXVm7JaFJIdJS+bTPHNbF0w46ETdjTXOZ/77huI7dD3BnbT2fpXMIwY2wTAPNWbGbb7o7ESN+jOhZlzvjBPBbw8vvxXpogs9icCUMS+Zf865yDY4L43numArln0Weaoe/xleP2TgQ7eHjp0IfVV7HPiAYezfFyqipd8WANDJygDo/qgOSfkLw3/mSn6XgZj9OzOExurk+suQKpdn2A3501ixs/PZe5k4ay3+hGHnwttfMJKiOToPTypUGww3pwXSUzxzVlbYf3W7/AP7DFEeQHT3Sugz+jQVXaQG1ofRUzxgzikTe6l+FPPZMr20VY5k4ayrPLNnUzj8bTRmGTm+uJRYXX1rZ2qwskk836OWTSEDq6lOdXZA4cyaSxQvI6eaax215ITkQ9efpIpo1q5KHXst+LQmCFi0HS4+7TO/aKaITzjpjM0Cz5xI6dOoKOLuW+V98J3B+U82vupKGowhOuxuMJl3T7NDiRRROG1nL3K2u77ctVjp9ZLU3UV1Xw6BsbWLe9vdtkUIDjp41g8brWrNFWXVleGnB8EQBVse6P7kxXE3w6xwg5lxbWMrSW+792ZCIgAVKXuz1p+kieW7E566S9Pz+1Agge7YOjNXlJFIfUVgYeM2PMIEY2VnPvouB7D76Mx2k3JhqRhGYM3TWK5oYqDpvi+EbeO2M0L6zcmnFejT83WRDfOXlq4nOm+R9H7zOcl1Ztyzpj3Z9bzKvjMxccy1eOdQIEPnBgclmmIG3vuKkjWPD21m4+ymzzX3rKCdNGsHNPF0+maeLpGkVlRYS9hjew4G3nme9IG3AGvU+HTBpKXWU06zuZyecCsPeIBqpjEZ5bvhmAqOufrKuMcuVHZnPc1OHMeyv781sIrHAxiBeTf/bBLfz2QzMDTQa5mDN+MOOG1PB3XxoMP0GRPAe2NDG4Nsb9r65DVfnZ3U60SnVF9/JFhNNnjeHpZZuy5v/SQB0pSSwa4bC9hnHzc2/z2trtjBjUfSb0B+eMo7G6gr88vSJjObk6Mi9YICg78ayWwYweVM2t897uts9PItdWllFsNCLc+9Uj3PM2pSx9e+rM0ajCHS+tzvRzfvAfJx9WtnHyLee9i4+/azzTRjcG7o9EhJOmj+Th1zdk1Maymay8awVkzQb8XndO1T/mBz9jniknk+YypK6So9zotvOOmBR4jBfGnek5huD5TSMHVSeeBRHhGyfszezxg1O0JY8zDxpHVITrn1wRWP8PzRnX7Tc95dDJw2ioruD2Bamz3YNCw0+YNoJnlm3m67e+lJhg7HFwmi8JXFPjfiO5e+E7ObMvBwmnyooIR+09nPtefYd4PGmaHe76eN9/4FjiSs73pLdY4WIQb17HtFENnD6rZ4thRiLCmbPH8dSbmwIjoYL88BXRCMdOHcEDr65LyTOWqcN+/ywnQsobbQeRS3MBJ2GhR1CajfqqCk6dOZq7F76TcR2ZbLZkgOOmjWDfkQ186KDuHUU0InxwzjieWLox6zyBXALMY1h9Fct/egr/+sK7U7ZPbq7n4IlDuP7JFTnn7mRbQ33a6EZ+eNr0FMGVzkcOaWFPV5ybfWuZ+EmaR7q35cNzWzh8yjB+fPr0lIikdMYNqeXYfYdzw9MrAu38XSGE8fWfOIjlPz2FvYY3BO4fP7TOHXyszLiMb0eG6Do//3vMFG77/KHB2QIaqzl91hhueX5lijD2HPqeplYIKisi/M/ssdz1ytqUCaLe6+h/ts50n9XbXliVyCB96dmzeOUHJzAnQLiA8y5t293Bv14MHsBkyi3mcfL+I1m3vZ1H39iQuN7e/KoJw5x7ceMzbxldEtkKF4N46ni6Iz9fPnboBBqqK/jNfW9kPCb9GTv74BZa2zu5ypesMhMtQ2s59YDR3PD0WxkXksrV6QMc7XOqD8tg6vvEoRPpjCuXPrgkcH+yIwsuo76qgnu+cgSnzQwW1h9913hqY1F+fV/mCXWecMk0EveTyel//tF7sXZbG7c8H9zpe+SanZ2LvYY3cMy+w7nm8WWBjv2EeSTgejVWx/jLp+fy0UPG58wW/fmjJrNlVwfXPNZ9Xk027cjDW4clG587cjJrt7UlQrDT6UybHNsTPn/UZNo741zzeDInmWdB6O17mM4nD51IXJXLH05m0QhaEnpMUw0fdAdenqm6JhbtFinm59DJQ9l/zCCufOTNQO0lm1kMnOjMliG1/OLexYnMEp/3+fm+eMxerG9t5/qnss+j6g1WuBhkT4Ee6kE1Mc4/ei8efH09/305VQ3vymBzP7CliVktTQkTQa5+9MvHTqErrlx4+8JAp+ydrv0328gyGhGe+s4xvO+A0Zw0fWTgMXsNr+fsg8dx47Mrmbdic7f92UbiYRhWX8VnjpjE3Qvf4QE3QubhxevZ0d6JqpOy3Ut305v060dMGca79xrKL+9Z3M2c6A++yDW3IwwXnLwvO/d08RM39byfTPc/X+ZMGMJ7Z4zi8keWsiQtnDeMcAnDYVOGcfiUYfz2/je6+XficSdKMJfmkovJzfWcMXMM1z6+PBH6XKhBXjotQ2s5Z+54bnzmLV5xQ7kzLU/wo9Om09xQxa/vdwaImSIEPUSEb5+0Lys370oRXh7ZHPrgtPXrJ+zNa2u3s2jNdma1NCX8leD4ZY+bOoLfPbAkZzaAnmKFi0E8dbQy2vvwx3MPm8gB45q44LZXEilJdrR3JuLnK9KG+iKSiJ4C+Pf5qaaddCY11/PV4/fm7oXvcF2azRqS+bRyjfZHN9Vw2dmzuoVV+/nWSfsypqmGL978YqJjfnHlFg77+UOJbLM9FS7gjNCmjmrk3Bvm8YErn+KT1z/PF256gUseWMK+37+Ha91Z773pkEWEn54xgy5VPvuX+SnmpN8+kNQwj5s6osdleEwZ0cD5R+/FP19czY3PvJWyL31NkN5w4fum0Vhdwbk3zEvRkgolwAD+74z9AWfJAX/iVm/+RzYTYVgufN80htZX8oWbXmB9a1tikFdVgHOn8/UT9mZEYzVf+Ot8tu7ak9DI029HTWWUrx63t+977rocNmUY7z9wDJc/vJSHX0+N7sqUW8zP+2aMToTVB+Vd+78zplNbGU34ZAuNFS4GKZTmAo7wuPKcA6mrquCj1z7Liyu3cO6fn+e9lzmTBoM6/dnjh3D9Jw/irIPGsfeIYFu4n/OOmMTJ00fy4/++yrVPLA/UYPJZHTcTjdUxrjjnQHa0dfLha55h+cad/PGJ5azaspvjL3kM6N0ouaoiyhXnHAgk84A99saGbqa4dIGcLy1Da7ns7FksWrONj137XCJnk9/P4mVw7i1fPnYKR+/TzPdvX8hffALGW2K4EB3/8IZqrvrobNZua+PMq55mzdbdtHd2cdEdXlbl3j/H44bUcsmHZvLq2u188vrnEr63XEED+dBUW8mVH5nNhtZ2zrrqGV51J7wWWnPxyrr8nANZt62d91/xFPPc5y0oeOesg8Zx+JRhiJCSEy0bPzl9OlNHNfL5m+anhIuHyQ8WiQiXnj0TgNluWLef4Y3V/PHjB/GrD84IVZd8scLFIIkRU4Ee6tFNNdx47lxqKqOcdfUzPLMsaVbKFrr7sw/MCBWpFo0Il3xoJidMG8GP//sq5/1lPvPf2swvfLmR0tP395TpYwbxp08dzLbdHZx62RMJE5ZHUNh0PkwcVsez3z2WfTII1ce/dXSvzu9x7NQRXHb2gby8ahun/O5xHntjAw0FmAWeTjQiXPmR2Ry1dzPf//dCvvq3BfzlmbcSmY17a07ymD1+CDd86mDWbWvjlEsf55L7kwK5tpcZhT2OmzaCS8+axUurtnHq75/k9w8tYb+L7gV6L/A9DmwZzJ8+eRDb2zoTK1mmz40qFAe2DOaPH5/DMt9SDUGaeyQi/PHjc3jkG0eFXrSstrKCP3/qYKYMb+AzN8zjF/e8zrZdHSnLpWdjr+ENLLjweH58+vTA/U7kXWYrQ2+QXGk/BgJz5szRefPmFfy83/j7S/xj/ipuOe+QxIS4QrBpRzvf+sfLPOhTla8450BO2X9UQc4fjyt/fGIZv7n/jZSkhoNqYrz4/eMLNhENnDVYvnTzi7zgprs497CJLF7XysWn7x+4jGy+qCpbd3XwzLJNfP6mF/jdWTM59YDRoZZDzoeX3t7K125dkJLp4K+fmdstz1Zv6Yorlz20hMseWpqSkPGli07o9frwfpZt2MFX/rYgkRZm5rimnKbVfHll1Ta+ePMLKalwLj5jOufMHV+wMlZv3c03bn2Jl1dt5YlvH8PgLOba3vLa2u3c+fJaTp81OmPUXE/Z2d7JD/+zKGVlToAVP3tPQcvJFxGZr6pzAvdZ4WJGuMTjyqTv3gXAP79waGK2caFQVXa0d1JbWcGiNduYNqqxYKM+j9Vbd/PHx5fx0OvrmTaqkcvOnlXwMsC5Vne8tIbaSie+v1xp6+jipmdX8qenlnPy9FF895SpuX/UQ5au38Gv71vMxh3tnDFrLB+eG5yfrDd0dsX594I1LFnfyrdP3LeggwqPto4u7l64lk079rC+tZ3PHzm54ALACeSI91obLgWeWrqRT1z/PHu64vzrC4cyq8D9Sr5Y4ZIDE8Jl8849HPhjZ3Gs/37xsJTUFRaLxdIfyCZc+q3PRUROEpHFIrJURL5T7PL90RmFyMRqsVgs5US/FC4iEgUuB04GpgFni8i0YtYhRbj0A3XcYrFY8qG/DqkPBpaq6jIAEbkFOA3oPgutF2xv6+DVNdvpiiudcSUeV7bs2sP9r65LWevBai4Wi2Wg0V97vTGAPyvbKmBuhmN7zOJ3Wjnr6meyHnP4lGFWuFgslgFHf+31gsJaUiIXROQ84DyAlpaeRdrsPaKBv35mLhWRCNGIk9o6KkJFVIiIMKy+Mms6fYvFYumv9Ffhsgrwp80dC6Qk5VLVq4GrwYkW60khg2piBZ/HYLFYLP2BfunQB54HpojIRBGpBM4C7ujjOlksFsuAoV9qLqraKSL/C9wLRIHrVHVRH1fLYrFYBgz9UrgAqOpdwF19XQ+LxWIZiPRXs5jFYrFY+hArXCwWi8VScKxwsVgsFkvBscLFYrFYLAXHZkUGRGQDEG71nWCGARsLVJ2+pL+0A2xbSpH+0g6wbfEYr6rNQTuscCkAIjIvU9rpcqK/tANsW0qR/tIOsG0JgzWLWSwWi6XgWOFisVgsloJjhUthuLqvK1Ag+ks7wLalFOkv7QDblpxYn4vFYrFYCo7VXCwWi8VScKxwsVgsFkvBscIlD0QkaBGysqO/tANARPrFMywi/TaJbDnTz96VoralX7yYJhGR/UTkKAAtYweViEwVkXdBebcDQET2F5GvA6hqvK/r0xtE5F0icg1wUF/XpbeIyEwR+YyIjOzruvSG/vLOQ9++99ahnwF3RPx74BhgJfAscLuqzhORSLl0aiIyCPgVcDCwAacd16vq0j6tWC8QkTuAE4ETVfUREYmqaldf1ytfROQzwJeAK4DrgY4ybUcM512ZA7wGtANXq+qzfVqxPOkv7zyUxntvNZfMDAYagKnAOcAm4OsiUl9ODxnwLZxBxAHAZ4GhwIQ+rVEP8ZmOHgN+B/wEQFW7ytQ81gL8P1W9UlXbylGwuOwPDFLV2ar6EZx+pRxTozQB9ZT/Ow/wTfr4vS/HF9IYInKgiOztfh0EHArUquoG4DZgM3C+e2zJ2mLdduzjfr0SuBBAVd/EeYH276Oq5Y3blimQWGE0gqO1XAOsF5Fz3X3xUr4nkNoWd2Q5HXhORI4RkXtF5Lsi8n53fzm0xXtXuoAzRWSQW/9DgGNFZJZ7bMm2xV0Kvdr9OoQyfech0ZYa9+s19PF7b4ULiZtyJ3A58BcROV5VlwFPAV9xD1sL/BOYJSKjS9EWm9aOG0TkGFVdpaprfKP+3cCbfVfLcKS15UYROQYSPpYFwNs4mss3ReTvIjK2FO8JBLblOFXdBmwHbgROxzGNrQUuFJEDyqQt3rvyEvALnDb8Afg/YBzwIxHZuxTbIiITRORu4I/ATSIyzTUZPQZ8zT2s5N956NaWG0VkH1V9y33vK93Div7eD1jhkjYK+QawQFXfBdwOfMrdfh3wbhGZqKqdwDqgDaihRMjSjn8D5wb8ZAxOx1xykVZh2iIitcAoYCKO6WIEMFxVV4lItLg1zkyOtnzG3f594ABgjarerqrX4yzNfVox65qLkM/YBTj+lv9R1b8AvwWWA+8uXk2zE9COZ1X1WOBh4IciMg34E3CIiEwq1XcecrblxyKyn7vPM7UW/b0vqc6lyFRD4ibtBDrc7Y3AayKyF/Ak8DyOYwxVXQiMx3FYlgqZ2jEIpx37QMKkNAXYrKovisjnge+LSFMf1DkTudoyVVV3AZ3Aczj28WOAFhGZUWI+i2xtWei2ZSXOaPODvt8Nx9GYS4lcbZnmjurbgQ8BqOomnA7t1eJXNyNeOzwtfhGAqv4ex/F9FrAG553/hbuvFN95yN2WD4vIcNcfuRd98N4POOEiIseLyP3AL0XkTPeleAKYIiIvAicBUeCvwJHAT4FRIvJ7EVmIs+7Ltr62v+bRjhtF5ET3Z+OBg0XkYeBU4BZV3doH1U8hZFsqgOtF5BQcW/gsVf2sqr6AY1ve2kfVTyGPtvzFNVteACwXkZ+JyDM4dv9FfdYAH3k8Y38WkROAe4ATReRXIvI4jhBa1lf19whoRyeOL2WWiBwgIgcAC3G04SiOuXWsiFxWSu885NWW8ThOfIBJwEFFf+9VdcD8AXvhhOSdBszCESDfcPftA/zTd+z3gd+7n0fgOPpO7es29LAdl7qfz8F5EI/r6zb0sC0XAb/2fY8Akb5uQwGer0ZgX+CEvm5DD9tyIXCJ+3kmTnTSGX3dhgztuBn4Ak4k6PeB/+IIzDluG7/i/q6k3vketuVL7u/65L3v8wtWhBuS6IDci3yFb9+ncEa8I4BmnPDWqe6+w4B/4M4F6uu/ArSjlDrhfnFPbFtK8xnL0Y5Pu+1odr9P8u07HzjX/VwS96VAbYn2Rd37tVlMRD4JrAJ+7G56BThbRCa432M4avuPgVYck8SXROTLwFXAA0WtcAYK0I4HAS0Rtb5f3BOwbaF7W/r8GQvRjgqcqKlL3O/L3d+dh9NZvwClMTO/gG3pG19kX0tmgxK/Hiea5cvuRd7X3f5bHHXySZww0P2Bu4E6nMlTXwT+DBzS123oT+2wbbFtKbF23AmMcPd/BceJf1Bft6E/taXPL6LhG9Ti/v8Z8Df3cxRn1HWY+32c+4JU9nV9+3s7bFv6vs79vS15tONPQJX7vbav690f29KvzWLqhHqCI+0nisiJ6qiI21T1CXff53DCK0spjDWF/tIOsG0pVfpLW/JohxfSjjrh7SVH2belr6VbEUcBnwUe9X0/GGfC5F3AyL6u30Brh21L6f71l7b0l3aUa1sGRFZkcTOaisg/cFI6tOM4IJeok3enLOgv7QDbllKlv7Slv7QDyrct/dos5uHemFqc2c9nAytV9Z5SvjFB9Jd2gG1LqdJf2tJf2gHl25aBtPrdF3CiLo5X1VJL5ZAP/aUdYNtSqvSXtvSXdkAZtmVAmMUgqVr2dT16S39pB9i2lCr9pS39pR1Qnm0ZMMLFYrFYLMVjQPhcLBaLxVJcrHCxWCwWS8GxwsVisVgsBccKF4vFYrEUHCtcLJY+QESaROQL7ufR7gQ5i6XfYKPFLJY+wE2b/l9Vnd7XdbFYTDCQJlFaLKXEz4DJIrIAWIKz8NZ0EfkEcDpO9tvpwK+BSuCjOGk/TlHVzSIyGbgcZ+GuXcBnVPX1YjfCYsmENYtZLH3Dd4A3VXUm8M20fdOBD+MkJ7wY2KWqs4CngY+5x1wNfFFVZwPfAK4oRqUtlrBYzcViKT0eVtVWoFVEtgH/cbe/AswQkXqc9d3/7lv4sar41bRYMmOFi8VSevhzR8V93+M472wE2OpqPRZLSWLNYhZL39AKNPTkh6q6HVguIh8EEIcDClk5i6W3WOFisfQBqroJeFJEFgK/7MEpzgE+LSIvAYuA0wpZP4ult9hQZIvFYrEUHKu5WCwWi6XgWOFisVgsloJjhYvFYrFYCo4VLhaLxWIpOFa4WCwWi6XgWOFisVgsloJjhYvFYrFYCo4VLhaLxWIpOP8fmAOCP0lBQFQAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Start a model instance, in this case a GR4JCN model emulator.\n", + "model = GR4JCN()\n", + "\n", + "# Model configuration parameters that will be used for the simulation for both the reference and future periods.\n", + "common_run_parameters = dict(\n", + " area=all_properties['area'],\n", + " elevation=all_properties['elevation'],\n", + " latitude=all_properties['latitude'],\n", + " longitude=all_properties['longitude'],\n", + " params=parameters,\n", + " rain_snow_fraction=\"RAINSNOW_DINGMAN\", \n", + ")\n", + "\n", + "forcing = 'reference_dataset.nc'\n", + "\n", + "# Run the hydrological model using our forcing data from the reference period, and with the common parameters that will also be used for the future period simulation.\n", + "model(\n", + " forcing,\n", + " start_date=reference_start_day+dt.timedelta(days=1),\n", + " end_date=reference_end_day-dt.timedelta(days=1),\n", + " **common_run_parameters,\n", + " )\n", + "\n", + "# Now plot the results\n", + "model.hydrograph.q_sim.plot()\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "7356dc3c-76b8-41a7-9210-9cad4a35ca21", + "metadata": {}, + "source": [ + "### Now redo the same thing but for the future period." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "4d6e838f-f311-4a6b-a3c1-d7301be3dd17", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAEiCAYAAAA1YZ/LAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABxNUlEQVR4nO2dd7wcVdn4v8/u3l5y03sICYEQQm+C0kFAVNDXAnZFUMT2wwY2bLzyvpZXUVERRVREUVEQpIn0HjqBQBJSSO/JTbl1n98fM7N7dnZ2dnbvzt69N+f7+dzP3Z2dcs7MnPOcp5zniKpisVgsFkslSQx2ASwWi8Uy/LDCxWKxWCwVxwoXi8VisVQcK1wsFovFUnGscLFYLBZLxbHCxWKxWCwVxwoXi8VisVQcK1wssSMiS0Xk5Aqf870icmclz7m7IiIfEpEHB7schRCR6SKiIrJdRM6PeMw3RWSHe1wq7jJa8rHCxTIkUdXrVPWNg12O3RkRaRCR34jINhFZIyIX+X4/SESeFJGd7v+DjN/misgdIrJBRKLO5O5Q1avc4+tF5K/uwEVF5HhzR1W9FNhvQBW0DAgrXCwWS7l8A5gF7AGcAHxRRE4Dp/MHbgL+AIwErgVucrcD9AI3AOcO4PoPAu8D1gzgHJaYsMLFUi0OF5EXRWSziFwjIo0AIjJSRG4RkfXub7eIyBTvINdk86qIdIrIEhF5r7H9QWM/FZGPi8hC9zw/ExEJK5B3DhH5vnvMEhE53fj9wyLyknvtV0XkY8Zvx4vIChH5ooisE5HVInKWiLxJRF4RkU0i8mVj/4SIXCwii0Vko4jcICKjKnNrM9fYS0TuE5GtrkbwZ3e7Z1ZKGfveKyIfzT1cfuIeu0BETopwyQ8A31bVzar6EvAr4EPub8cDKeBHqtqtqlcAApwIoKovq+qvgfnl1FVVe1T1R6r6INBfzjks8WKFi6VavBc4FZgJ7A181d2eAK7BGf1OA3YBPwUQkRbgCuB0VW0DjgaeCbnGm4HDgQOBd7nXK8aRwMvAGOB/gV8bQmmde8524MPA/4nIIcaxE4BGYDLwdZzO9X3AocAxwNdFZIa776eBs4DjgEnAZuBnQQUSkWkisiXk7z0F6vJt4E4cTWEK8JMI9Tfvw6vufbgUuDFM+InISLcezxqbnyVritoPeE5zkxc+hzVV7TZY4WKpFj9V1ddUdRNwGXAOgKpuVNW/qepOVe10fzvOOC4NzBWRJlVdraphI93LVXWLqi4H7gEOilCuZar6K1XtxzHdTATGu2W7VVUXq8N9OB33McaxvcBlqtoL/AmnY/6xqna65ZwPHODu+zHgK6q6QlW7cUxK7whyNqvqclXtCPn7Y4G69OII6Umq2uWO6qOyDkfL6FXVP+MI3DNC9m91/281tm0F2ozft5KL+btlmGOFi6VavGZ8XoYz6kVEmkXklyKyTES2AfcDHSKSVNUdwLuBjwOrReRWEZkdcg3T9r6TbAcYRuYYVd3pfmx1y3a6iDzqmri2AG/CESAeG12hBI7GBbDW+H2XUYY9gL972gfwEo45Z3yEMkblizimp8dFZL6IfKSEY1f6tIzMMyrAdvd/u7GtHeg0fm8nF/N3yzDHChdLtZhqfJ4GrHI/fw7YBzhSVduBY93tAqCqd6jqKTgaxQIc01PsiEgD8Dfg+8B4Ve0A/uWVqwxewzHvmRpIo6quDLj2NHHCbgv9vTfoAqq6RlXPU9VJOJrSlSKyF7DD3aXZ2H2C7/DJPh+V+YyCrrUZWI1jgvQ4kKwPZT5wgO+cB1Cmj8Uy9LDCxVItLhSRKa4d/8vAn93tbTgj/C3ub5d6B4jIeBF5q+t76cYZDVfLeVsPNADrgT7X0T+Q0OdfAJeJyB4AIjJWRM4M2tE1i7WG/F0XdJyIvNMIhtgMKNCvquuBlcD7RCTpajQzfYePAz4tInUi8k5gXxxhGsbvgK+KE5QxGzgP+K372704z+rT4oQsf9Ld/h+3rCJOUEe9+73RFeiRcc/b6H6td89RrvC3VBgrXCzV4o84PotX3b/vuNt/BDQBG4BHgduNYxI4ms0qYBOOL+YT1Sis6//5NE647GbgPcDNAzjlj93j7xSRTpy6HjnQcvo4HHhMRLa71/qMqi5xfzsP+AKwEcep/rDv2Mdwwoo34Pi93qGqG4tc71JgMY4J7T7ge6p6OzjRXDgBDB8AtgAfAc5yt4NjJtxFVpPZhePnKYWX3eMmA3e4n/co8RyWmBC7EqXFYqllXG3vZaAL+IKqFjWNisilwEU42meL4RuzVAkrXCwWi8VScaxZzDKsEZFfFHCK/2KwyzYUCAkqOKb40ZbdGau5WCwWi6XiWM3FYrFYLBXHpqIGxowZo9OnTx/sYlgsFsuQ4sknn9ygqmODfrPCBZg+fTrz5s0b7GJYLBbLkEJElhX6zZrFLBaLxVJxrHCxWCwWS8WxwsVisVgsFccKF4vFYrFUHCtcLBaLxVJxrHCxWCwWS8WxwsVisVgsFccKF0soW3f1csA37uCJpZsGuygWi2UIYYWLJZRnXtvCtq4+rrh74WAXxWKxDCGscLGEYpf1s1gs5VAzwkVEfiMi60TkhYDfPi8iKiJjjG2XiMgiEXlZRE41th8qIs+7v11hlz2tDDZ5tsViKYWaES44a2+f5t8oIlOBU4DlxrY5wNk4y7WeBlwpIkn3558D5+Ms2Tor6JyW6FjRbLFYyqFmhIuq3o+zTrqf/wO+CJhj5zOBP6lqt7tG+CLgCBGZCLSr6iPqLFTzO5x1vC0DRLGqi8ViiU7NCJcgROStwEpVfdb302TgNeP7CnfbZPezf7ulTMR6XSwWSxnUbMp9EWkGvgK8MejngG0asj3o/OfjmM+YNm1amaXcfbA+F4vFUgq1rLnMBPYEnhWRpcAU4CkRmYCjkUw19p0CrHK3TwnYnoeqXqWqh6nqYWPHBq51YyHrc7HCxWKxlELNChdVfV5Vx6nqdFWdjiM4DlHVNcDNwNki0iAie+I47h9X1dVAp4i8zo0S+wBw02DVwWKxWHZXaka4iMj1wCPAPiKyQkTOLbSvqs4HbgBeBG4HLlTVfvfnC4CrcZz8i4HbYi34MMezM1qHvsViKYWa8bmo6jlFfp/u+34ZcFnAfvOAuRUt3O6M9edbLJYyqBnNxVLbWJ+LxWIpBStcLKHYUGSLxVIOVrhYImEVF0tc3Dl/DXfMXzPYxRgyPLlsEz+7Z9FgF6MoNeNzsdQmNv2LJW7O//2TACy9/IxBLsnQ4L9+/ggAF56w1yCXJByruViiYVUXi8VSAla4WEKxiovFYikHK1wskbDzXCwWSylY4WIJxS6HY7FYysEKF0sk7DwXi8VSCla4WEKxiovFYikHK1wskbCKi8ViKQUrXCyhWMXFYrGUgxUulkiodbpYLJYSsMLFEor1uVgslnKwwsUSCau3WCy1Ra1bE6xwsRTBqi4WSy1S47LFChdLNGr9RbZYdjfSNd4orXCxhGJ9LhZLbVLbosUKF0tEav1Ftlh2N2pccakd4SIivxGRdSLygrHteyKyQESeE5G/i0iH8dslIrJIRF4WkVON7YeKyPPub1eITY41IOzNs1hqk1pPJlszwgX4LXCab9tdwFxVPQB4BbgEQETmAGcD+7nHXCkiSfeYnwPnA7PcP/85LeUQ8zDpkcUb2furt7FlZ0+s17FYhgtWc4mIqt4PbPJtu1NV+9yvjwJT3M9nAn9S1W5VXQIsAo4QkYlAu6o+ok6c3u+As6pSgWFKtRS/n9+3mJ6+NE+/tqUq17NYLPFSM8IlAh8BbnM/TwZeM35b4W6b7H72b7cMkLgHSQ0p51Xs7u2P+UoWy/DARotVABH5CtAHXOdtCthNQ7YHnfN8EZknIvPWr19fmYIOQ7wbGvd7nHQ1pHRttxeLpWaocdlS+8JFRD4IvBl4r2anpK4Aphq7TQFWudunBGzPQ1WvUtXDVPWwsWPHVr7glpLwrG+13mAslkJ09fbT05eu2vVqvamEChcRmS0it4nIrSIyU0R+KyJbRORxEdk37sKJyGnAl4C3qupO46ebgbNFpEFE9sRx3D+uqquBThF5nRsl9gHgprjLOZzRzP94X+WMcKn5JmOxBPPOXzzC+3/9WNWuN9TTv1wFXAn8AfgPcDswEvg28NNKFkRErgceAfYRkRUicq57jTbgLhF5RkR+AaCq84EbgBfdMl2oqp6x/gLgahwn/2KyfhrLAIj7PRbXAFfj7cUyBEmnlXf/8hHuWbAu1us8v3Irjy3ZVHzHClHrTSVV5Pc2Vf0ngIh8W1X/5G7/p4h8s5IFUdVzAjb/OmT/y4DLArbPA+ZWsGi7NVUbHWU0F4ulsnT19fPYkk08u2ILC759+mAXp2Jo9SxwZVFMc0kan3/o+62+wmWx1DDxay7edax4sVSWhAxPrbjWTcjFhMvPRKQVQFWv9DaKyF7Av+MsmKU2UN//uLCJFCxxM+yES43XJ1S4qOovgV4AEWkwti9S1c/GWzRLLRG3RlGtkOdqcPUDr7Jo3fbBLobFZbgGi9R6baKEIv9CRBpxHPuW3YyquVyGieLS15/mO7e+xNuufGiwi2Jx8d7h4TBwMal1E3KxUOTjgHnAA8CTInJsVUpl2e3IaC41Px4Lxyt9Z1df6H6W6jO036x8an3CcdRJlMNkXGkpHecNjt2hP0ydrpbaodZH+qVS6wOxYj6X+4DDgWOAw9zkkpYa4KFFG5h+8a0s27ijKteLfRKld53abi9FGerlH84Mu0dT4xWKorl8XFV3AZ+IuzCW6PzzWSerzYOLNsR6narZq4fJPJdaH03ujgxbn8tgF6AIUYTLz12H/s/iLowlOnVJ59H1V8nwWq0XebiZLiyWuKj1plLMoX8s1qFfk1Qr0WO13t9M+pcqXS8uar3BW4YPta4lF9NcxPffUiNUe0Z77PNcsuFiFktFqfVOuFyGdLSYdejXLpnoqpivk7FXx3wdO3qxWEqj1k3I1qE/xKna6CX2UGTvMrXdYIpR4+19t2S4PpNar1exrMioape7ZsqnRGS6eYyqvjXGsllCyPpc4n3DqjU6Gi4p94e6cLTUPiJDo50UFS4u/8BJf/9PoMYTPe8eSJUNSfEnrqzOdSy7H8P1nap1ARNVuHSp6hWxlsRSEtWOFquWQ7/WG0wxhnr5LbWP4LTLdI2/bFHTv/xYRC4VkaNE5BDvL9aSWUKpdi6u+K/iBSjEe6UXVm5lv6/fzrrOrljOX9vNffek1h3f5VLrtYqquewPvB84kaxZTN3vlkGg2iP9+HOLVec6v314KTt6+rl3wXredfjUeC8WI1feu4hZ49o4Zc74wS6KxYeqxro+kbhOl1oXmlE1l7cBM1T1OFU9wf2rqGARkd+IyDoRecHYNkpE7hKRhe7/kcZvl4jIIhF5WURONbYfKiLPu79dIcN0FarqhyIPDw2pIeW88t19/bGcv1oN/n9vf5nzfjevKteqFnHdu2p3wXFHcA6VKWFRhcuzQEeM5QD4LXCab9vFwN2qOgu42/2OiMwBzgb2c4+5UkS8JZl/DpwPzHL//OccFlQ70WO1ljmOGxs4ULvU+EA8MnGnZBoq/smoZrHxwAIReQLo9jZWMhRZVe93Q51NzgSOdz9fC9wLfMnd/idV7QaWiMgi4AgRWQq0q+ojACLyO+As4LZKlbNmqNK8kGppLFIlaRl3lF2Nt/eaZrjcu7gd7eK69GvdLBZVuFwaaykKM15VVwOo6moRGedunww8auy3wt3W6372bx92VHteSPyai1Ofak0Kjas+1XgeTy7bHP9FBoG0KskYhH+1++Dhlky2XCIJFzcNTC0R9AZqyPb8E4icj2M+Y9q0aZUrWZWomidpmC1zXK3Jp3Gys2d4rnI5hB9JDv1VsiHX+v0qlhX5lmIniLLPAFgrIhPd60wE1rnbVwBmqM8UYJW7fUrA9jxU9SpVPUxVDxs7dmzFCx43wy5xZbWvE9cFarzB1zKxmZOq/EzScftc3P+1ng2imObyBhG5OeR3AeZUsDx+bgY+CFzu/r/J2P5HEfkhMAnHcf+4qvaLSKeIvA54DPgA8JMYyzdoDLuU+1WKfoubajT4Ce2NsV9jMKj1kXhU+qxDHyguXM6McI6eShRERK7Hcd6PEZEVOH6ey4EbRORcYDnwTgBVnS8iNwAvAn3AharqxZZegBN51oTjyB9+znyyPoof3PUKnzppVuzXq95iYfGePyPEarxhhjGEix5KXIK52iP8uDWXzHVq/CUOFS7V9LWo6jkFfjqpwP6XAZcFbJ8HzK1g0WqSavkoqrVEbLVDhOO6TjXae433KWVT6+uTRCVun8tQSfIadZ6Lpcao9tzQuEd/2QZTLZPC0J2wV+u29nKJ7ZkMs2ixoTIt3AqXIUr1gsWchlItzSVuqp1NOg5qfcRaLsNFc0nHnDe+2hOoy6Vk4SIiI0XkgDgKY4lOtUcv1VqJcqjnMKtG9F6tdyplM0zqFXsoskuta7CRhIuI3Csi7SIyCicVzDVupJZlkKjWCLxaHVn1NBeH+JzH8VPrnUq5xOWgrvbd6uzqjfX8nkm81jW9qJrLCFXdBrwduEZVDwVOjq9YlmJUXXOpVhRX3L6dIRLGGcZQLnsYtR79FJW3/vShWP0u1Z7jVi5RhUvKncT4LiDOSZOWiFR/gn61JlHGepnYGerlH0zii+Cr/kPpi9PxMkSSr0YVLt8C7gAWqeoTIjIDWBhfsSzFqLbmsmF7D+f9bl7sMfyx+3ZinqxZDZPVcBVgw0VzAejrt+9BVOHyT1U9QFU/AaCqr6rqf8VYLksRqhWKbI767npxLdvisidXyVw1HDSk4epzedvPHo7lvINxt+KcpZ9t+bX9HkTNivyCiKwFHgDuBx5S1a3xFctSq3T3xaPuZ+a5xN1g4l6qwE6iLJuVW3YNdhEqRl9/fGaxoZJlIpLmoqp7AecAzwNvBp4VkWdiLJelCFWboe/73htTo/Hqs76zO3xHS42PV2uPweiEY9Vc3LYyLKLFRGQK8HrgGOBgYD7w5xjLZSnCYE0GjMuW7NXmmoeWxnL+7HXiHfVVJRS5ipmwl23cUZVreVRrLZS4qYZZbLhEiy0HPgvcpqpHqeoZqvrd+IplKcZgrecSVxRMteqTiPk61Wjv1ep/r3/8NY773r08uWxTbNfwd5A9MZldq02cZjGP2hYt0YXLwcDvgPeIyCMi8js3U7FlkBisJCY9fXFpLtWpUcakMKRHyNUp+1PLNwOweF182otfGHf19gfvOJBrDEI33BtjtNhQ8blEXYnyWRFZDCzGMY29DzgW+HWMZbOEUD2fS+4bPNQ1Fw8bilwb+KsSV8BItYl1notLrUcNRhIuIjIPaAAeBh4EjlXVZXEWzBLOYPlc4hyRefSnlWRM9ivvvsWWaqQa0WLxX6Jq+M1icWgulbhhNz2zknFtjRw1c3Sk/Xtj0vBh6ITTRw1FPl1V18daEktJVC05nt/nEle0mPG5tz9NMpGM5Toe1TCLTb/4Vr56xr589JgZFT1v1fK9ederojjr6otBuFSAz/zpGQCWXn5GpP17Yg1Fdv7XunCJ6nPpEZEfisg89+8HIjIi1pJZQhnRVDco141NczHsYnFG2ngdZVzC2X/Wax9ZWvlrVKlXqYap0qvJIdM6AOjqrXynPBh9cFwh+w5VmhM2QKIKl98AnTi5xd4FbAOuiatQltrB34/1xuVzMT7HGWnjya245Je/449DDlS7S4lTlnnnbqp3NNVYzGKDgE3/El24zFTVS920L6+q6jeByur6IYjI/xOR+SLygohcLyKNIjJKRO4SkYXu/5HG/peIyCIReVlETq1WOauJ+WI9v6J6yRJim+diSJc4/Tqer2UoR4vVeqdSCt7ou6kuPuEyGPcrTs2l2kuCl0tU4bJLRN7gfRGR1wNVydUgIpOBTwOHqepcIAmcDVwM3K2qs4C73e+IyBz39/2A04ArRSReA/4g85afPhjbuas1Qz+VMM1i8TVMr6OJa7JedRz6Rq7qIS5pvOI3uMIl7mixag0qtnf3xXZur6XUeqLPqA79C4BrXT+LAJuAD/l3EpG3RzhXl6r+K3IJHVJAk4j0As3AKuAS4Hj392uBe4EvAWcCf1LVbmCJiCwCjgAeKfGaNc1g2VvjS/9iCJcYNRevMx6sgIjKnDT7sT+tpJLxRtZV4061uGaxHTF0ymZb+fHdC/l/p+xd8Wv4+dT1T3Pi7HG0NETtYqMj2UiLmibqPJdngANFpN39vq3Arr8CbiJ8jt+xQGThoqorReT7OFkCdgF3quqdIjJeVVe7+6wWkXHuIZOBR41TrHC3WcrAPzKuhi05TpPCAws3AEPbtGQWvS+tpGLSy6s592hCeyMAq7d2xXqdhxZtqIpwAejs6otFuCTcBxNn4EslCK25iFxUYDsAqupf6vg2Vf1IkXP+oZQCur6UM4E9gS3AX0TkfWGHBGzLewoicj5wPsC0adNKKVJNMFidYzUmh8XZaBas6QTgtw8v5dw37MnUUc0VPX9VzGLGNdZ3dle8DmHXi+vcTfUpRrXUs3zjzvguRnUH+3ENktob61i9tSv25ZQHSjGfS1uRvxxUNazTj7yPj5OBJaq6XlV7gRuBo4G17uqYuP/XufuvAKYax0/BMaP5y3GVqh6mqoeNHTu2xCLtPvgbY88Q11xM/v3S2oqf02+ujMMnYl5j3RDPIu3VRQSmj26OJe2++Qji9FO0+bSUuCLfvGkIW3bWtnApprM1q+qXROSdqvqXgVxIRE5R1bvKOHQ58DoRacYxi50EzAN2AB8ELnf/3+TufzPwRxH5ITAJmAU8PpCy1yLVVlz++2378+W/Px9bmLDZCVfD9AaQSkaNZymfOGpiKnbbdsXXwWSjkuL0gbnXAkY217NmW7xmsXJkS9TBTtLn+4pjzg5Aa6PTbW+N8dlXgmKt600iUofjPB8oZeUhU9XHgL8CT+GsJ5MArsIRKqeIyELgFPc7qjofuAF4EbgduFBVh0fw/CDgNcZ9JrQC8TUYk2qY3gCSMTgV/J1XHG4LUxCv2BynGSn+BIneqUVgRHNdLKNxLfA5Kr+4d3Hkfd//uj0yn3fFpLl4gZVbdvbEcv5KUUxzuR3YALSIiOnEF0BVtd3cWURuLnAeAaIl5QlAVS8FLvVt7sbRYoL2vwy4rNzrDQmq7HRprEtSn0ywZVc8L7RZnWoIMIAqKC6xYD75r900n/cfNT2W61RjPoUnKAVhVHM9m3bE3GGW0W6ialOquUs6xD0hdNNQNoup6heAL4jITap6ZoTzeRmTt/u2C044sGXIkW38Hc11bNkR/wsdewfjkkpUXrr4uy6JI+SqtoOEykIExrQ1sKu3nx3dlY2yMjW9cm5dKY/QfN5xaS5edTbtqG1/W9TWtcC/QUT+J2C/R4Gdqnqf7+9e4OUBlNPio9r9i4hjE98ckypu1qdawiWOuS7VmNRY6zmlSsGsyZjWBgDufyW+HLnlPJ6oGci9Z3/jJ44G4tdcNm6vbbNYVOFySsC20/0bVPV0Vb0n6ASqemwpBbPUBmZj7IjJJu5nY4zCZd+JWUtuNWZrx6K4VDkrcpwXNE89prUegAuueyq2a1RDMI935+zEJVy8GlRrEFYuocJFRC4QkeeB2SLynPG3BHguygVE5M2VKKglF6/BvGGvMYxuqY/9erFrLm592htTsar7E9obGNfmjJBj0Vx83+MULsfMGkNdTLPzoUo5rLxoMRHGus8l1suVUZmoSwt5p/aeSdxh+5t39tR0+p9ihs0/ArcB38XN3eXSqapRF9b+FnBLGWWzRKCjuS7WVCbmmUe21LF5Wbyay8iW+li1IyWbxywOxaWai4UlExLrhFOpSrSY59NzgkbipiyzWAkjBJGsLy8uzdgTKL39Smd3H+2Ng7P8RjFCNRdV3aqqS3Fydqnx1yoiUae1D9Zy77sFyYRUx7yD0NFcH1v4o9fJ1CcTMecWy85HqNZ9qzRe55JKJFAd2skrM/NcBOpiCLDIu14VTp6MmJ6lq7efb9/y4oDyqW2t4YixqCEZt+LcOgEacVKxvIyTebgYHyuvaJYwvA4lKRLb2iTOdbKfm+qS9KWVdFpJxLQMcdyjcUdzcTqxeDIj554zzmAxTwOLK3lldsXDaghhqEvFPw4tpy6lHCNI5MHL7x9Zxq8fXEJTXZLPn7pP9PIYn2s5v1ikoYKq7q+qB7j/Z+GEFQfmeReR77j/v+UeO+xmx9cSiYTElj7eRITMuvZxvNBe+00lhf5YU+5rph5xpAKpZm4xrxOLq4OpRvJd89x1MU08ynHol1GZqIdkzJURNRfv94EswFeNtl8uZT1NVX0KOLzAz0+IyM9wUrRYYsIcvcbrc8meO2mMlOMilUjEPhrzGn+tr4dRGKfcdTEKSTAT1MZyevfcmrlWXMJloES9v6qaMwgrdly5qfPN09aycIlkFvNlR04AhwB5wegicikwCjgH6BORg1T1W5UoqCWYREKqZrbImGFijLJKxayJqZpCMobz+77Hk/7F+e/lRqvlDqYYXslFHH9bPNfQwM9RKeX2CoaGX8R3WIlFv2r52Ud9mmYm5AYcH0zejH13+WOAo9zvVrDERMY0IvF3xh7eOhL9MTrc4/e5ZP0T1TCLlTtD//YXVjP94lvZuD0/LNu7PV7Ia9wdTKxmMc+hD7GGVXuUE2AR9TXxdvPckcUGYYkKaIa1LFyiLhb2TQARaXO+qj+9i8lvVPVlEflNJQpoCSeZcBz6jkoe75yHVMbGH8eQ342ASgrdMeYWy9Vc4m+YqTIDH37z4FIAFq3bzujW3Pkf3ug77npUw6Gf0SREMvWp+DWM4pfTREpZAkLEGVAkE8V9h15ZSn185u7VWlG1HCJpLiIyV0SeBl4A5ovIkyIyt8Du73T//1clCmgJxpzrAPHM2fCT0Vxis/FDMmafi2q8Phe/2aXcDtMT4EHHZ8xiXtRbXM+jirMIhJjysFUAb5XMYpiPwREu4ftnfFoD0A1rWXOJaha7CrhIVfdQ1T2Az7nbgrAO/SoS9+jVPGsqxmtVzeeCZoRkHPNc/P18ucLFszwGChf3f5zPA0zNJZbTuyfP/frGOeMjd+blXKIcAdbRXOceW3xf7/yOubqI5uKVr2SHfvaAWhYuUee5tJg5w1T1XhFp8e9kHfrVIzPPJeaIoSwS2VFZ/hWq4HNxZ2slY46y8yhbuIRqLu6zj9nnUg09wnTog7PKYpwKTDmn9l6TYq+LqYGkImguiYzw3r01l1dF5GsiMt39+yqwxL+TdehXn7jDas0XvxrzQ1IRbNUDug6uEJPijb+s81fIoe8J8LC+I5PGJublb6qzEqU74o9Zcy3ncZTyvnunT0TyuQzcpD0chMtHgLE469ffCIwBPlxg32uAA4C/AYjIV0XkRhE5ZIBltQSQqJJzOu5JlM41JHbNBXXqkkjE66j+/Bv3BqC+zAgo73kGLSvt97nEtXJnNVwgmdxi7rUSCan44MV8zuUJl6jXyX6OMv8so7kMxOdSww79qNFim4FPR9z3GRH5naoeICJvAE4Fvg/8HDiy7JJaAsloLlVYvHEgCfl6+tL845mVvOOQKYGpY7wGVg2fi5AgEVMIt1ePvca18Ya9xrCzp7y8UV4HGyRo/dFicZtE451E6fz33oi4Q+vLoaT7awjJovUoU3PxIh770xqrlj9Q4poS6y1kcAbwc1W9CYg/L/xuhJkuBeIfwTj+EOdzOZrFdY8t44t/fY7rn1he5BrxJ64UcTuxGM17jnZUft4377Cge5H37GPqXzyzzeNLNsWm5fl9LnGYxcyzxfpuGZ+jDJKyPpfSrxXnROBKEZdwWSkivwTeBfxLRBoGci0R6RCRv4rIAhF5SUSOEpFRInKXiCx0/4809r9ERBaJyMsicmoF6lOzJGL3uWQ/JweQ8HGzm711fWfwWi25Ppd4OwARxxkeZ1Zkwek8yn4u7mFBJi/vlHUxm8U87l6wjlufXx3LuTPpXwyfS5yKSznvVinviVePhBQ372aXNCitTIoakYK1K13iEi7vAu4ATlPVLTgRZF8YwPl+DNyuqrOBA4GXcNaXudtNpHm3+x0RmQOcjZOx+TTgShGJf6GIQcLTJuJOH+/4Q5zP5TRQb5QWdqjX6ccbLaYIQipu3w4DM/F4QilIOHlbkjE79E3j5bKNO2O5RqZ6A9Bc1mzt4uoHXi3YSZuby3nm5iGh7cz0uUQYvAx3zSXU5yIiPyEk+4OqBvphVHUnjuPf+74aKGvoIyLtwLHAh9xz9QA9InImcLy727XAvTjrzpwJ/ElVu4ElIrIIJ4vzI+Vcv1bx7O5xT2zMTVzpjZTLES7hozRva1WixYT4fC4Zs5jEZhbziN0kWsU5jZkoqzLMlRdc9yRPL9/CSfuOZ88xeTMkMtQlpSwtzxTwfWmlvkB4uaJZ814UzSUz4CpRc9FspGDcWutAKKa5zAOexFnD5RBgoft3EFm/StzMwEmSeY2IPC0iV7tzbMa7QssTXuPc/ScDrxnHr3C3DUuqlcpEGNikvUSEhiTEHy1mmt/iXcXRNYuVeQ2vnEH32j/HKS5hbM7Qj83nYghjcDTxUt+vXT1OV1R4zXovWKQ8f55Z96ipYJIRot6yM/RLxxvo1XJm72IrUV6rqtcCs4ATVPUnqvoT4CQcAVMNUjiC7eeqejCwg9wll/0EDSvynoCInC8i80Rk3vr1eQmeax6zk4T4TCOBiSvL6DCLxfRXIhHnk8s2cfZVj4R2AIpTj7h8Lv4lCspt/OHRYg6pmE0jZpLiuPowc5ljKO/5expcsY4/VabJ1TwkTDipGvVISGxZkU3NpZbNYlF9LpNwMiJ7tLrbqsEKYIWqPuZ+/yuOsFkrIhMB3P/rjP2nGsdPAVb5T6qqV6nqYap62NixY2MrfNx4Yb1xj2DMxJXlCRfnf6FiKk7LHIjP5Qt/eY5HX90U6h/w1tyIa90YM1pMBhCR5h0W+Fw9QVxmgMXHfj+POV+/veh+SWPZ4biUPPN+gfE+l3DBZJGBi0cqIYHzhophPoOeIsebUW/F2mQi0yhKLlLsWmsliCpcLgeeFpHfishvgaeA/46tVAaqugZ4TUS8dUBPAl4EbgY+6G77IHCT+/lm4GwRaRCRPXG0rmG3Gmbe6LUq0WLl23mL+Vwga3qLewkBL81M3DPBkzJw7SgwFNlbLKwMYb9oXSd3zF/Lzp4IVm3jWcU1Sz/jzzd8FVDa+5zVisMd+qlkeQMK85BiWrFHFPNuuT4XiD8MvRJEnUR5jYjcRnYS5MVup18tPgVcJyL1wKs42QESwA0ici6wHDcbs6rOF5EbcARQH3ChqlbLP1R14kzCaCLIgFLNFD3WGI33p3VASwiEHeY4Xb1oscq3TH8HU7ZD371PwT6X7PmhtI74tU27opch4JqVJi8U2RCYdRFjPIvl6PK21iWE7QP0uRQ3dUVPY1Ou4qLokNBcoq5EKcDJwAxV/ZaITBORI1S1KhqBqj4DHBbw00kF9r8MuCzOMg02A+lgSrqO8XkgiSu9hhQ20nLMVdnOJVVm6pSwW+FpLvFFi2U7S5HyAy28w4Keq19rLWlgUcItNS8d19ClkOZSygAmEdUslkyU9TzMI8LMYjl5+CK8X9l5LiUXKavh1Vg2A5OoZrErcZJRnuN+7wR+FkuJLCUR91wHDzO3WHnRYkVMF+7/AQnLCB2nN0O/XOduKWVJSvkOfc8MFejQzwRzlB4aXoq4zil7FcyuUN47Vkx7NzMalBWKbJy32PGm7yiq5nLzs3ku4VC89C9OeYa+cDlSVS8EuiCTa8ymcxlEqpVfyhyNDcS/E2VVQ29yIwxsRBZ2DU+viMvnYp4xMRDh4mkuISPlcjriUkyN5lnj68PcJyLZme1Q2mAp6oqOdYkEaS3dhJzjc+kr9m45RPEdDmRxtDiX6q4UUYVLrzvLXQFEZCxQu8a+3YhqqseV0FwKHemfuzGQEVkx4eeZ32KdRIk3ei3zPO7/KIkr56/aGvm8ZndWbO5KrlksXs3FDOGF0gYwxbViZ3t2me7S6mLWvWi0mPs/ynpB5YoWZWATmqtFVOFyBfB3YJyIXAY8CHw3tlJZijLUfC6RJlGaPpcBJBgMExqmEIunYWZH4skBpPX36hCY/sUw8wD85D+LIp/XHCz/7J7w48xONa4lj/0+l3KWkPAipouN4lPJ8nKx5aSPCfW5ZD9H0Yz9v/71yRX8/ekVkcpUlr+tykSNFrtORJ7EcaALcJaqvhRrySyRSFbpJfNCa6Fcs1jESZTJgY/IijtSHX/FroIzussnR3MZwDyXbrds4ZMoo6cGPPNnD3Hw1A5OmTM+s+37d77CJ0+cVfgg0xwUU8xr9n65UVZlaOLescXerboytWJTaPUWG/TkLHNcTDPM/f3zf3kWgHtfXs+Pzz445MChkVss0tspIr9X1QWq+jNV/amqviQiv4+7cJbiZCdRxnQB47ypAWRFLj6JcuApZryxdVjn4Tn0EzFpLuZIfCARad5RgRqcz4QYhWdf28JvH15a0mJZaVWa6pK0NaaKmoPKxb9YWCY5aknzXNxzFdVcytO8tQwhOxCf3k3PFHfwZ4Vw7UqXqEOf/cwvrv/l0MoXx1Iq1fK5OIkYnc/ldMpR0ot7K1E61yi/0UQJAU0lYkr/YozEkwkZcJBVmOZSV0aodinmLU8Q1yUTVdBcHMqZt+W9MwUHLp7mUqZZLFdzKcHnUlRzCblm2ADJXFhvqDr03XVROoEDRGSbiHS639eRnRFvGURijxYzVJeBrETpnafoLOoBaC6L1+8AiuR/cjPXxuVzycxzEcfPFEvKfZ+/rRRKOcTTJuuSEholNRD86V/KSTEU2aFfps8wd4Z+ofc3d3s0n0vh36OkmYlL+64UxRJXfldV24DvqWq7qra5f6NV9ZIqldESQH5m3Jg1F7Ja0kBSaITJQKEy0WJhgtYUYnGYFDJmMQa2HrxX/eB5Ll5nWfpyTP7w19CwbXX2j1Vzyblj5S0hEWWtIMg69EtpK129/fxnwdrM90L3wS8ko0SLhf0cmmbGeIeHg0P/Enelx1k46fe97ffHVTBLNGKPFjMjYDKjytI7Gs2MxAv8Tm6HWaqwzEnREWpSyPpD4kxciQxsnovXaYTdh0poq1t39dLRHDxlLa2OllefTMTncwnolKHE1R8j5hari5g92eTy2xawdlt29dTiZjHXoR8hK3K4cCli2pVoa8YMJlHTv3wU+AxOhuFngNfhLL51Ymwls4TivZhREkJWgtzULKUf7xUvtJxSfriz2QkXE36ePySeSZSavUYZDn1VZfmmnZnBQmBuMfd/OT4X//3v6i18rx5ZvJHOrj4mjWiKTXPxyPgqBqK5FLnX5QxcVm7JzcVW0Czm+x4lO0PYrz19xRNkJmpcc4mqV38GOBxYpqonAAfjLOBlGWTiDkk0X93EACJUwnwIMHCfizmCC19zQzNCrJz5OkUxRuJ1SWdGeCl1+ePjyznue/cWWSzM+T+mtYGmuiSTO5oin99/ujCh8fLaTgDqUlI8BLdMspqLaxYbQPqXYgKplEmUyzfu5Knlm9mwvTtne1HNxfAdFbtO2EArioZUjaW6B0IkzQXoUtUuEUFEGlR1gZEC3zIIeK9UJdKlREGQzEi5p4yOJuNzCb1G1vRWakRPruYSbhYD114dg7ZnehDqUlkzTDIRLcXv08u35HwPug/Zawgnzh7HgjXbSihfbp2jdE7V8Ln4NZdSHn/URezqXM0lynIDx37vnsDthX0uuddORFhuIVRziXC/B+LTqwZRNZcVItIB/AO4S0RuImABLkv1GcjqkFEw310RcezvISp74fOE+1w8yg2tNjvJ0BGseoEDcS8W5twrKM3G7zd0hS1zjJSegNN/a6KUrSqhyD6fSznzXIppxV6HvWpL9GUH/BT1hbj/I4UJGz/3pzXzvkARs1gmoCPeNYkGSlSH/tvcj98QkXuAEUDxpewsseG9t/Wp0juwcvAacEOqPOFSzCzmXCObuLLUjt9MyxElYWBsucWMSYF1GeFSuonHI8x0l5mDUsLz8N//KO9NfTKebAZQmfQvkhmQhO+355gWIPtcyqHQu+8vbSIhRVMYmVpkb3+aRAJmjGzh1Q07Ipnf4gpKqRShwkVERgVsft793wpsqniJLCXR4AqXcjr8KPjNKPWpBN19pXc0xRz6Aw2t7o/qcyEb8lzOkrfFMCcF1pWhufiji4OOzb2G0FvCvfLvGsXvVJcUtnXFpblkAyCgPM21mEPfe4czbaXM597WkGLTju7A3/Lm65QYitzbn6avX2lpSGW+FzzO/Z9K1rZDv5jm8iTZ9uhHgRkVL5ElEv4G012laJ7yNRf3f4FDlew6K1CG5hLV52IIsThT7pt1Ke1+5Ta1oI4wqx0JqUSiJCHp11yi+LbqyjSFRiFTGr9ZLAaHfv0AB2JN9Um6C2ou2WcC0SY4mr/29KXpSytN9Y5vrtB1TIZ0KLKq7lmtgljKY6ANphj+9upoLuWbxcJmJXu+ECg9C0CO5hJhnktckTZZzaw8n4t/Bn3YcxVcn0sJZje/5tgTMvN+VEs9p+43nm1dfVVL/1JOxonsKqfhPpeGlNNxl9tWGuoKv/v+4kaZ4Gges2zTTgCaXeESZko1MzTUskM/6jyXY4O2V3MSpZvPbB6wUlXf7Jrs/gxMB5YC73IXMUNELgHOBfqBT6vqHdUqZ7WoVIOJTMbnkizTLBbu0PeHIpfa8Zsho1HW0aiG5lKOz8WfWDKoMzOrV5dM0FtCaJV/12KaSyqRoD6ZiC0UGd+IP5O4sqLpXxyyJuTy/EdR3n1zuebimkv297df+TBgCJci7dnLwxffcxk4UUORv2B8bgSOwDGZVXMS5WeAl4B29/vFwN2qermIXOx+/5KIzAHOxkm2OQn4t4jsrarxeCQHmdg1F9/3hrpE6MS7QkRK/2Ikrix1Ls3/3L4g8zls9UYzpUlfWkmnNeNErgg+fwiUGi3mM4sFCRdvX08DK0Vz8X0POzatSkJcv06VE1eWM4mymEDKBr+U1yGHmYTzlg4wNPBC71dQcZvqovtcGlLxRfFVgkhhE6r6FuPvFGAusLbYcZVCRKYAZwBXG5vPBK51P18LnGVs/5OqdqvqEmARjjAcVngvWDLhdMg9/fHKTq/RNNUl6SohckhV+d4dC3hptTMXI3JywZLX3Mh+DjeLOb/VD9C5W+z8IkJdGZF8eWaxUIe+kHKFZNQMDf77H1b/dFqrkFvMwRvxZzIXlxFhV9gspjnnLveZN4SYhAstHRCqVQY8M09zibLiZX0qQXcZA71qUW5M3gocAVMtfgR8kdyllcer6moA9/84d/tk4DVjvxXutmFLuXNPIuFrAE31pQmX7r40P7tnMbe9sAYo7kMoN1rshH3GZT6HO/Sd6ww0cijs/OAIifqyzGJRNJdsR1af0Y6iXcN7dl978xygWDYDp+OO1aHvG/GXE/1YLLeYR32qnACLLA2pZNHO3Ht63uClWOSin+aGCGZu9YJ5krHlfKsEURcL+4mIXOH+/RR4AHg23qJlrv1mYJ2qPhn1kIBtec9RRM4XkXkiMm/9+qGdyaa+zAiuUvD6vKa6ZElzHvwdfSHBlPW5lD5yBRjX1pD5XFTrESPKrsIjv7TRWZYTiuwnULgY1St16d4f370QgEkjGose55nF6lPx+VzMJQog2ymX4tfLmsUKXMP972hh5Zn4nvjKya5DP/z99fB8oWHBL0GysDljFiviNxRPk6pda39UzWUejo/lSZyElV9S1ffFVqpcXg+8VUSWAn8CThSRPwBrRWQigPt/nbv/CmCqcfwUArIJqOpVqnqYqh42duzYOMsfD8abuXVXL397amU8l/F9X75pJ6+s3R75eL/tvCukMYiYmZdLNYtl9w9daEmdjj82s5jRWWbT5ZRiFiseiuyREzQQcb2VDZ1O4EMU/0NanZDauqSwq7efzTt6Il2jFLKxdQ5ROmU/xRz6Ht7co3IGYmPbGqhPhpnF3GuUICSDTJlN9cWDDjI+l7ryzGI7uvu48t5FsftrovpcrjX+rlPVh2ItVe61L1HVKao6HcdR/x9XsN0MfNDd7YNkFy+7GThbRBpEZE+cZQIer1Z5q4nZD23v7ovlGv7sywtdwRLZxu/r6HcVyOvk7VW+zyW7f5TOPK5ACLPU2Y6/FBNP9nMyIaHHCpKx0e/sjfb8S8ns4KXcX+8KpI//IarxIDrqky5ZjbJ0zaWQtpsXXVeuz6UuZJ6LbzJoFM040CxWH1FzobylEP7w6DL2u/QO/vf2l7nxqRUlHVsqUc1ibxaRp0Vkk7EiZfRsefFwOXCKiCwETnG/o6rzgRuAF3FS1Fw4HCPFyjFSdPX2s66zq6RjPA3Ca8CfOXkWEN3Gn28WC2sM5S9zbLaxMJ+Quh1mXCHc5kztskKRjc/NdcE2dVOwe8JlR3e0V9wTwvUZx3m42SYh2aWal2zYEekapWAuUQDOaBxK01w8n0uYVuztV59KlJV4FVwzVCGzbuYa3r7FHfOBZrEIDn1zGkKpmstX//FC5nOhgV6liGoW+xGOdjDaWJGyvcgxFUdV71XVN7ufN6rqSao6y/2/ydjvMlWdqar7qOpt1S5ntSg1gPbca5/giMvuLukYrzPywim9EVlUv0ueWayIz6VcjcK8TrGyeaM+KM22H41sZ1lWKLLxUAuNlE0B1uqmC9nZE01zSfvvc4RQZO+YWJIkGnUB87mUnnVgZ0HtPVvugQS/NKQSrNralRP2nrlCns+lPM2lqS7aoEdEQn1AUYh7cn9U4fIa8ILGvSKVJTLmk/jkCXuRkOKmqocWbSz5Ohnh4rb+Rvfljxox5ldACnf8jkbhNa5SR1We+a2pLsmunvAGLVK6EFu9dVcmnDqMIM2lXL9OU31wR2hakjwzSqmai3dc2HN0hItktMhYMhq4/z2ZmkomSCWkpE7Tq9Omnb3B18gEWZQ+Z2dUSz2vfOd0IPtu/vzexQEXca/htpNyfS7eO/Pgog1Fy1aqoPRfL+7Z/VEnUX4R+JeI3AdkpkKr6g9jKZUlEt6L3JByFqXq7ddMuGUYjmlI2N7dh6rS1lhXcF+vHSZ9wiWqOu43b+3q7c9cP68+OI0rlZCSs/B6DaWlITxU2nPolxryetR3/wPA0svPCN3PNI+Uo4WZ/V6jG2rqv19mWv+WBs8sVlhzMTsV7z411SURCRfiaXWuEedcCvV1ylB69KMn8/757Cp+cs7BBffznklp846ywR/e/yD869KUm1DWy0f31PLNRa8Vlo4mCL95Nm5VIarmchmwE2d2fpvxZxkkzNQRXoe/YvPOSMd6L9kh37qL/b9xZ+i+WbOYdy3nQzH7duZ437uvGj4xEEoPd4asyaa5PhXp2MzIMqZ5LoKUrOVBbmaCpvokqvkag9mRZTSXELOYebj37BMJaKlPFVw4SzMaa/ZZB4wHBox/8iGET1YMPEex7MPG57qQiK9ipEIyOfizIjfUZaPeVJVPXPckd85fE3iMSV0ywSHTOjh65ujQsjhztZL0pTWyudIvVEvJglAOUYXLKFV9u6peqqrf9P5iLZmlKN6r7nX4J/7gvkjHeS9ZFHNNOuPQdzWXVGkdZtAL3FXAbOU1zMYSJ2qCqbmkwjUX1/wWX7RYtrP0nLOlmPhMQeId76+P2ZF5mkvY6opB5o+ECE31SXYViDLzipEQ4aJT9gbggCkdkepQCv70L1Cao7qrt58XV2XNlWFaiSC0NKRKjKzM3ruwNEH5IdVZs1ha4V/Pr+H83+dG2wW1jVRC3InKURz6pb3D/kFKKVGM5RBVuPxbRN4Ya0ksJWG+l94oKSqlmAW899FvFouaXyxoVBWkWeRpLiX6XLwqtTaEaz1ObjGz8ZfWwIrN9fCv55JKCDtLEJSmIGhtcMyVfsGRNb2J4XMJM4vlb/OEXyFfTdrQXA6Y0sG+E9tzVkqsFP4oK6CoozqdVq59eCk7e/q44A9P8uyKrZnfir1b7Y0pOruiC5ecCauhmos3qPD7XNIFfRte27jpwtczsrkuc1yU998cIEX1T/nb4jOvbYl0XLlEfVsuBG4XkV01FIq822OuDlkKpTiY+zU7EodsqGiYdvD4kk1cee8iILezbCxyrJm/rHyfSyq0YXoGpXJDkQ/+9l2hv/s7y6b60gSlOVejrbG44MjMcylHc6lLFjwu7ess45oNbi5R4FHMLHbni2u59Ob5fO+Ol7n3ldzsGl0h90EE2hrr2N4d7PgPwtRWkv6V3AzyQ5GzwqWQ2cqzCsyZ1J55H+tTCRrqkqFmZ39G9KgDJL//8+4F6wrsWRmiTqJsU9WEqjYNZiiyJYv5ujaWrLlEt7WqG44qPrNYd1+aLTt7eGBhfuqcd/3yEf739peB3NHS/77jQKDA6NL0IdUn2VWiE9nrpFrqw81iHg0ljvqism6bN4/IuV/N9cnIYcKQe7+yYca+MhrCoi6ZoD6VCPW5BOGY1FIFy+afPFuqHyQqgZpLqvBkRYBtXY5wWLRue55WFqRRq+be01I0F1NZCVPc/OY9s+MvqLm425NGNvCGiJoLRlBKVBNikJCLMwA4VLiIyGz3/yFBf7GVyhIJb6RfsnApKXpJc1KSmNrH9+98mff/+vGCEzNVc52N3ki8kFaSzV+WCB2BBuFdprk+mlnMW/FvZ8QQ3qh859aXgGxdmkOc5kGYdvhCmosXTu3R2pAKrUchzaWlIcWOIppLIqOxJiNreZ1dvfz2oSXROq4An0uxpbS9d+qBhfnhumHPXnDuaWdXX+RO1Xz3wzWXXI9+ZjJob39BzaU/7fj/EgkxlkdO0FKfDDdzuv+zE06jvV9BGQziGDB4FNNcLnL//yDg7/uxlWoY8Jd5r/Gtf74Y2/nNtlGqLfyX9y+O/EJ6+aU8zAio+19xGnehUVZvf1a4/N+7DwwNBjDrEzXiK7ec0cxioJGjrAaCV59S/Uem07W1MVhz8TI7m/z+0WUFrxPUtwmOf6pQJ2Y69KE0zeW7ty3gG/98kXtfLp4Q1r88sHetQj69Hd19/ODOVwqeL1grztLWWEd/WiO/X6Zw2Wd8foDsRTc8wz+fzaYu9PZucd+v7d19BScr9qc148s8aGoH4AyOPKEfJgDNgJGog5egeUpxztIvtszx+e7/E2IrwTDlC399DoCvv2VOfBfxRjvJ0mJEr3/8Nca1NWa+96c1o5b78WZpe3iqeFdfumgesK6+/sxIvKO5PqMtFPa5OJTnc3GPjRBpI+Lk7QrzOZhEmTyZXx4vNDraNTz6c3wujpPXLwCdiLfsQ9nkBhlc8/ASPnH8XnnnDOqkUskEzfWpgsLFGxSYPoSoAxJP69wYIdFlULRYc32SLQUmRP70nkU5q476CeosTbeOpw12dvVlBhhhmBriG2aNYcaYlpz5Ljc+tZIbn1rJ418+Kee4ZEJoqU+yvauvYCLVfs0uJHbZ2/bnPUdOY1x7Iy0NKfrTSndfOtAqYZqAIdwnl3O9gJRKO3v7GRnp6NIpZhY7XEQmGN8/ICI3uan3R8VUJksETB9FOdMPNu/MNvxikTlJc1SZmUTZn2kYhcwl3b3pTMNKioTmwfL7kEoZUb26fjt3zl9Dwk2H0tOfjuR3aWlIRgpLfctPHoxcFg81hF0p0WI5ZjHP59JdXHOBwhF8QX1bfTJBa0Nh4ZKd5+JcaVtXH6+u38HzRmRWIbw1SRavL5492z8/BJxnWOi5FDPNhT13QQzhEs2p75/bM2dSe3BKnoD9WxudehSaT2K2rRFNdRw9c4xzXENW6ylYLhwtvdh+JtXWXIrZU34J9ACIyLE4ySF/B2wFroqtVJbiGB2MqXVEtSWbTv2wkX6/5vpcWlwB0dnVl2kYhRq808k7vzWkErS7I/Gtu/IbtqNRuNFi9YmS5rmc+IP7WLCmk2RCaG9yrrGtQOehZH1VzfWpkHxUWfyNMsoI3hOkzfVJdpXp0Pf8W/maS/CExkJLQwdrLs7s/kLml6xZzPk/b6mTuu9XD7xarAqZTi8wTYq/bO5/c3nnlhChV0xLDzS5YmqDTtm2RXTq+5edLpQ9wJw86+EFDxTUXNIEWgyyg7DwMkYRQiamz+Vjx80ABle4JI2EkO8GrlLVv6nq14B8/dtSNcwO5qCpHZnFn6Laxc3srsXSpZg+l1QyQXtjii07e7KaS4HQ5u7e/ox5q7k+xYimwsLFpNSllD1EhHav89hVeETu3TdnQl3p1wlznp8yZzwAU0c1A6U79E1B5plf8hz6mt/pQchiWQF9WyohOeYXP55Zz+v8vLJEmaXfVEKASTbkObuttTFFZ4EOsy7EqQ5FHPoCY1qdReXWbI2WHdz/7NoaUoFaT1CmgdbGOjrDNBefyTlzDXcQFhbVJpLVEIsJoasfeJWr7l+cM3B53QwnA0Cp5udSKCpcRMQzTJ4E/Mf4LWpeMksMqGqmgxERLjh+JlB4xO7HjKMPE0hOtFjutlEt9Wza2ZsJzSxoFutLZ0Jdm+qTNNY5YbPbgjQXY3TZXO84NEvNwpsQImouDi0lhgl7I8WwIIBRzfVMaM/6s0qd55J2I4j2GtfK4dNHUZeU4I62BM3F68D3GteaPVwk1Gbvn+fiBY1EMcGObqnPfH5hZbgZzRtNm36MtoYUPX3pwPcqbJY8FOhojddotCtcipnFPA1934m5TvwRzfVs6+qjP605Gl+Q76itIcX2rt5wh35AfTqawwdh3rWymkv4+/WdW1/iv/+1IDNw+d1HjshYEQZTuFwP3CciNwG7cJY3RkT2wjGNWQYJ1dwY/CijHRNvTXsolrwwvwF0NNezZWdPJjTTr7l4+7/rl49kzt1cn0REGNFUV7DReKM+T0BEtYt7qJLRjoIEmJ/mkFDcIJpCfEYeaUMzAqeTKiUirS+d5vDpo/j3RcfR1lhHR3M9W33ObSW4tyoUWOFt/eBRezC5oynTkYeZVfzzXLzUL1EyI5uDgjcX8Vd52SLM2e8ZIR5QLr/AqfOZyTYXCAQAN0KuPuvQD6OjuZ4xrfV5iTA7jHfTvBVBPhcv7DnMoV+OcHHq4gSkJCS6Q7/PuNfZ7OPxREtCEeGiqpcBnwN+C7zBSLmfAD4VW6ksRfEy1np4tuTtEYWLqamHvcTpgAzGo1rq2bSjB69d+xu812A6u/q488W1QNaOXFC4GOXJCojSXvy+tIb6dSDXt1NsPoEfz+S2ZWfhKCglN3y1vbGOrt50SSk6zACKkc11OcEX3kXMJ3Lyvo4prpC7zdRC7v7ccTz7dSeT08gW514Fdcj+eS7/925nAuykjqbidShB4fSES50RTh/mqJ7U0ZjzvS6ZYHJHEzPHttDakMq/V+QGi2SzSIc/j7QqJ84eR0dzfc52r+PfsrM3Z8a7fyVKyAYmhM3Q9y9rDdDRVJ+5Rhie9hnV5+INcloaUiWHMZdD0QkSqvqoqv5dVXcY215R1adiK5WlKOrO1/DwRnthIzKvc/SzKSRkNJ0mp7MDp4Ft2dmbESJ+4WKOQv/jpphoKiZcyI76ovpm/PSnlfamcIetacpoaUgVFcZmxzC2zTGphN4vn+bS4ZqI/NpH2PVMp3VHc31e5+936F9xzkFA4UADMyKrsS6ZeRYj3Y5z04780F7/PJe2xjpGNtdFMiMWGqkH4QWWmMLFDBcuxOwJjrlqZ08/D37pBO76f8dl3stCiAipZILGukRRrdhby8ZPRrjs6s15NwKy2DCmrYH1nd0FO/9iZrEgQQm5mmtY8IMf796MaKrLWAdKbWOlUPlMdJaqoKbzANMsVvhlKZTgcmNA5+LRH+B0HNlcz+adPZnGV0hzMfFs9h1NwR2A2R2VK1yAjOZSyCxmdswjm+vYsit8Loa5aJPnDL7useUFk39682g8RjYX1g6C8Hc4I5vr8jQl098Gjo9qzzEtBUfjfhOXx+gWT1gGaC6+eS7edaKMdP2mszDfmTf6N81bXsLOIHOiJ4y8ztEpo5BISOa99OPX6EY01XHj0yu59KYX8vb18FsGssd6WkVPYL3MI/YZ30ZfWguGZBcSLo11SZrrk6GDGK9oa7Z18ZcnV0TKnvD3p1e6daijo6mOZELYuL34XKRyqXnhIiJTReQeEXlJROaLyGfc7aNE5C4RWej+H2kcc4mILBKRl0Xk1MErfbwkAsxihaJsoPCIMuwFSxsTvTxGtdSzs6c/04kU8rmYeA11dGt9QWHmdZgDES6NdUkaUolA09XDizbQ2dWXuU5Hcz1dvelQn5PZ6R07ayzgCJwr7l4YuL/6RrwdRmcUhb48s1iA5uITYOBlOC402z7XxJU5t2cWC+jEggRS1OAEf7qZMO3Q6xRThubiaZ9B5fL294S2yciW+tCs1aZmvGlHD9c+sqzgvhowqIJcf0iQ5mIKJG8wUqh9FfK5gNPGNhaYLBpk/owS/bZ8k7PeU0dzHYmEONcIGVgOlJoXLkAf8DlV3Rd4HXChiMwBLgbuVtVZwN3ud9zfzgb2A04DrhSR0pJvDQH85pf2CA79Qs7Y1zYVXmQsyC7sNbC1bpJG/6ipJWTm89i2BjZs78kTdKa5aiDCxbvG+s78RvOeqx/L+T7KNVkVMj845XL+n334VN5+yOTM9uUF7llaczvkjgFqLl7wRE5kEvlRWy31qYKBA+kAfwA4ptS6pLApoP7+ReIgerYBr9OdMaYFCH+OniZipjDyIrqCRu6e6W9EU4Bwaa4LvM/+AAjz2IL+kAABDlmH/uYduZrLa+5CfeYh3nXMetz2/Oqca/tNzh6jWxtCMxz4D4sSJbpso1tGb6DXUs+G3VlzUdXVnn9HVTuBl4DJwJnAte5u1wJnuZ/PBP6kqt2qugRYBBxR1UJXAf8sbc9RGWYWK9SQlmzcEbgdnEbmH12Ncm31KzbvAvLDGT3fRBBjWhvoT2teh26aqwYqXMa3N7IuQLh4ZM1ins+hcAPzRurvOXIaqWQiM2u+cKeU6wsb2VKa5uL3uYxsrqO3X/Oi2vwmm5aGwh1/0Cx47xwjm+vZFNDBZLWd7EEjmurYEvJMFq7t5OBv3cmqLc578aXTZwPhHV8mgsmosxfKHNS59rhphzynv8nI5mDNxR8mbAqXQj6kQj6Xkc31pBLCus7unHfgve7AJcff5g4szIHONQ8vzb1GAc1ltBs0E4T55v3h3COB8AGSyev3yq5wOaa1ITSVzkCpeeFiIiLTgYOBx4DxqroaHAEEjHN3mwy8Zhy2wt1WE9wxf42Rlr18/PmlUkknm2pYh+xfzwGc0eu6beE+F38b80fQbPB15H3pNEdMz2YHuuxtczOfPcETNmLy5sMU84cUYnx7Q0arCsKrThTNxROcXuimZ+svlNXA72z3BHHYKPSsnz2USX7oz0I9rt25X2Z9gswizSEpU4JMNpnyFTCNpAOOGdvWwPqQ+3rdY8vZvLOXW55bjUi0QUJQKHJjXZKW+mSgOamnL019KhGoHU8Y0Uhnd1/BAVZ28JJ9f8P8VEHCJZEQxrc3smZbV+DkSPOQ0a3OddYY92zOxOxKJWGay6gw4WL43Ma770dUzfjdh0/LKd9u7XPxEJFW4G/AZ1U1LJNg0NPKewtE5HwRmSci89avL569tVxMc0Z/WvnY75/kXb98pALnzbehj2tvDDQHmddvSCWYOCIbzjmpo5G127oK+mNU8xvAqJZc4bLeN/rp61c6musyDanVGGWOdU0e/nKampiIMLa1gfUhQi+McW2NoQLTGy2Oask3W/jxhIuXQPAn73HmPRTKSuD3uTTVJxnRVFfQJp5OK8+8toVPXf804Ahzs6OdOMIJ/V29xRAuvkhBcDScQvXIrBsS0Nonjmhk1Zb8smVzi2W3jWtrZP327oIphrykpjt7+kglJNKco960Up9M5Ak+xyyU/wx7+h3hcvYRU5k1rpX/fO64zG+T3TBpf338pTU160IC2W92NpkwopE1W7sCU9ibNNenaKpLFhzo9KcLTwod55p2A9dgISvEvIFeVM3YzJ4wtrWBdZ1dsa3pMiSEi4jU4QiW61T1RnfzWhGZ6P4+EfCWVVsBTDUOnwKswoeqXqWqh6nqYWPHjo2t7ObL4WkOSzcW9nFExTlt7os5ri18xN6XVj527AweueQkzjpoEuCsH9GX1kC7O+SPpCHfmervyHv709S5IZ+Q64PxGnbQGjBmBzNhRGPOiK8UxrU30NndV9Dk4dVnnDuTfnWIM9RLQ+OF7x4yzYkbuf+V9cE5udL55qeJIxoLXqMvxymsLNu4k9vnZye4ZjvMXcZ++deY3NHMlp29gU597x0MGolPHtnEqq278rb7Q5HBGSX39mvBUbInXDy/k2feCjO99PalA/OFjW4NHrn39KWpTyaYMrKZuy46jhljs1kHJo907tXKLYXal3OdCe1Z4VI4cWfw/QKnnT28eCP/Mvwn2Sv4hWR9Tps0hYAzQTm4pFNGNtOX1sD2bM7V8kxv34y4vIcpXCaPbKKrNx0pe3U51LxwEecu/hp4SVV/aPx0M/BB9/MHgZuM7WeLSIOI7AnMAh6vVnn9mAOPUtOZhJM/shrf3sjaAiP2dFpRzS549OmTZjF7QhvvOdJRkwsJJf96LpBvFsvTXFy/gefTMMN2J3U0IQKvbcrt0Px3pphpK4ivvdlZ3mC8u5xAoXvh0d5Yx4imutCABs9s4k06g6zmFpiAk3xh7AjK/A4cct+JL7rLNJgmt/HtjYiQIwCCQli9yYWmEPJfIxWQl6uQUAqKMPOWaSi0OJwZ6p5MSCZa6ms3zefxJZsCj+lLa84cF4/RLcFBGT196cwiWX6muILY8wV6+AcBI4zBUViEXaFMM54Q/e5tC/J+87fJ0a0NGWGckFzzVZhZbMrI4LqAZxZz8O5d1L6lqT5776aOdPLfhb3/A6HmhQvweuD9wIki8oz79yacDM2niMhC4BT3O6o6H7gBeBG4HbhQVeObhloEMywzSuqMqASlXZ8wwjFxBY2o/aaRGWNbuf2zxzJnkmO6Wh1gGgEvWix3W30qkdO5+X1Iff1KKpHgsyfvzaQRjRw6PRMlTmNdkgntjSzblBtEoD4zRJig9DNnYjtjWhv48NHTAZg22mk0SwsEKpjPZOqoJl4LaMAenV29OekyAL5zluNDCjInBc2PmDiiqeD9Nf1gf3lyBQDffOt+mW31qQRjWhtyjndCWHObrqfhrAwRLkGjZE8o+Y/z5xaDrH2/0HMxI76SCckZlPzp8eWBx/T0p/NSuABM7mhk5eZdee9yt6u5BDG2rYGmuiRLNgQ/d68qr3dT20O4WayQ5hK2IqUfU0sa09qQp7kUMotlhUt+x+/365110KSMwCuG6W/ykquGvf8DoeaFi6o+qKqiqgeo6kHu379UdaOqnqSqs9z/m4xjLlPVmaq6j6reNpjlNzuy/lLyYhQhSG0f19ZAd186cESd7WByH/meo51w0UINMii3GGQT+41vb2BbV19ORFBvvxPRs/+UETx8yUk5C5MBTBvVzHKfadDfmMe3N7K9uy9Saov+tHLoHh2ZhuqFwL66PrhOpq186shmVoSM3Lp60zTVJXM6Wc9ntTrAnGSOKs39N+7oCcz0HDTiNLUkcDoaM/S5vz/fnOKZhIIE3i3POVbhoM6yUCcWNM/Fu0ahka4pJPx1LZQZuLcvHai5TB3VTGd3X9673N2Xpj4VPLNARNhzTEvec/dfeVx7I588wUnq/sirGwMzGxSaRAnw4ddPz/nuZQuA/Hs8ze3AAfYc08KiddszEXJhmouXZsev4UN+299ngrPGTJi/1WNMa1a4TCnyPAdKzQuXoU6uz6VywiXI4eiNRJYF+HT6MqYRn/+kpZ4xrQ0sXNcZeJ3+Ao3M81fMnTQCgIVrs7OQPbNYIfYY3Zznd/KnmfFU9mUhYdLZ66VzTD6jWuoZ0VTHkg3BM6OXG1rT1FHNrNi8q2BAQ1dff54ZZlKAH8TDWaIgd9sEVxgFBRkEvRNNPuEya1wrC9f57q/vIuPaGkklJM/f0Nef5kp3XZWgQcKMMY7PYtG63HsVZBYb39ZIfSpRsDMy6+LNXzltvwmB++bUJeBdmZIx2eTeY8+hX4i9xrXyaoHnbl7l86fuw7i2Bq55aCmfuf6ZvH0LTaIEmDt5BO84dErm+1td/yXkty9TuBw5YzQ7evozz7I/XVhzaaxLMmlEY+A77A939wTGW39afFE7bz4cOKljxrTWR2pj5WCFS8yY0b+mFrNs446CKUSiEDSRbqbr3AxKN+FpTUEdjL/zyrmOKkFywutg95/iCBdzKeDu3v7A5Vk9ZoxtZcP27jwTgSnDZo5rcetS/MX3z8UREWaMzR/Beswy1kKfOqqZnv504LyYeUs38cfHltPgGymPbW2guT6Z1yFn6kFwB7M4oKMIijjyr4ey9/g2NmzvzoxMHXNK7jHJhDAhIPLruO/dm/kc1JGNbKlnXFsDC9bkDi6CggASCWHqyKaCWm5vQF2++/b9AViwOnjw0tNfSHNxR9U+jaqnrz/UBDRzbCsrNu/K1ZwKjOk8Z/jt89fwh0eX5Qww/JNh/UwyIi73Md4nf/uaagiXg6d2ANnBX39a84SRyeyJ7bwUcN/8qZ8O3cMxO4cFpnj434HPv3Ef3nLgpAJ7DwwrXGKmkM/luO/dyzf/OT/wmKO/ezdf+fvzoefVAI1ij9HNpBISKFw8237QKHHW+FYWrd0e7KspkLnVa1yjW+oZ1VLPcyu2ZH7r6ktnIsWC8LSdF1ZmBZK/MU8f3YIILC4g9EwczSW3jHuNbeXlNZ05dfJOf94xMzLbpo4M7sQA3vfrx9yy5d6XREKYO2kEzwesVdIfEACxn+vXeiFgieCguUd+4eIJw8Mv+zcL1mwL1FzAEfh+rcL0pRQywewzoY2XfcIlO0M/95j9J4/gyWWb+f4dL+c59oMGSyNb6jnzoEm8sq4zz1S2o7uPW59bHTgI8ASyX5D19KVDhcve41tRJUdYZhfyyq2LOYr/6j9e4BYj+ivMoQ+52aHNwYe/fZmayz6u+WzZxh3c+txq1m/vLpj+BRxz2+L12wPNduYAZsbYVg53/ZrXGpM0o3D2EdM4ZlY80bJWuMRMf4jP5f5XNvh357jv3cOqrV1c91iwA9TD7wAHJ3Jk2qhmFq/Lb6x96cKay9xJI+js7uPVDTvo609Hcjqe5KZ5nzt5BAdMGcFzbsfZ25+mP600FrCLg9HZrsp2tuobjTfWJZk6sjnSOuz9/fllPHBqBxt39ORE28wc28oZ+0/Mmaezx2jPP5N/Ha/TCBqR7z9lBPNXbcvYzz3SAaPRtsY6Zoxt4bkgYRSU/ND3YPc1bPqPLN4YGGQBTmf08prOgia+Qh3Z7AltLFy3PUc4eB/9dTls+ig27ujhp/cs4nM3PJvzm/9eeLxxzgRU4ZW1uQIsbAJfW2MdU0Y25WjE4JrFCsXvAoe4o/inlm3O+81fe3/U41p35K+qgYM3E299G3Du64yxznvkF/re+wWOH3Fkcx3zV23jwj8+xbKNO0O1o30nttOX1hyTs1c+/6P0BNefn3gNP22NKaaNaubGTxxd8FpxYIVLzORqLrmNz9+x9Kc1x18SZjbzR4x4jG6t5/b5a3LsqP1pzXT+fhMPkInmmrd0E9+65UUO+tZdmXxh6XT+ZE1wlvN9+muncPC0kRw4pYNX1nays6cvb9JhECNb6pkysimj7fz2oSXcvWBdXkObNa41z1zj58anVrBqa1deJ3iQa4J4+rUtmW19/em8DnaPUc20NaZ4NkCr8EbIQc/hgCkj6O5L55kt+tLpQA3hwCkdOdpddv98QeBF8HmMM1a23LarL8/H5LHfpHY6u/syzn//iLdQR3bItJH09KVzNDHvXfUfYzqv/RNDe/o10B9y4FRHU33GeBZQfMnk/Sa186IhXOav2soLK7eF+lzGtzcyuaOJecuyoc+F5gh+yI0u9PBm9hdKl2Oyj3EfUknh48c6K8FOGJEbvJJMCCfvO4437DWGZEKYNb6NJw3BF6a5eILynF89mrM9KO/ZF97opNoxVxvN7J9WTpkzPjNHq1pY4RIzhXwukD+r1p+2opBtG/LTrnt4a2P/7amVmW2/euBVzvvdPIBAc9WMMS2MaqnnkcUbud4NGfXmroRlbvXyZh04dQRpdcxcXRnhEv5qHbHnqMwo/BvuBDD/SPHgaR0sWrc9dC2Ui9zRs7+Msye00VSX5Al3fsWWnT0s3bgzJ4U+OGafA6d08MzyLXnn9gRkkHA5aqZzn+9fmJvdIZ3Od+iDY05au607b+6Of4Dx2w8fHpiU0esI//TEcu6YvzZw4uN+nrnR1Qi/c8tLodfyOGJPJ1XPo69uzKkH5Jt5TH+VP0lqr6tV3P7ZY/i7MUqe3NHEpBGNPOS798V8jnMnjWDJhh2ZdvLtW5z3pFh6+aNmjuahRRvz6uvvkN8waww/dTMuQDajuCfwg3xBHub7VpdM8M7DpvDwxSeyt3F/PK7+4OH8/lwnveHB0zpyTJVhWpgXXu5fzdLJ0JBbmRHNdRwxfVTwPKeQNhwnVrjETH8BnwuQl4zwaV8Hd9mtL3FDgJoL3oqH+dsvOmVvGusSrDTMQa8Yo/8gc5WIcOgeI/nHM6syJiBv7kpvAYeriWciePa1LXS7EwDDNBdw0tdv3tmbYxp71jeyzZg4luebOCC3k/FrLqlkgtfvNYb/LFiHqmac70Gzvg+a2sHLazvzhLsnIIM6s3FtjcyZ2M79r+QKl0JaxQFu4MP7f52bmdl06P/wXQdy/D7jCOIbb92PY2aNyThtg7Jf7z2+jca6REag+icuBgktcCb67TO+jYcXZYWL9976NRfzHP6EpX3unJXZE9o52Bgliwgn7Tue+1/ZkON3KZRCx+PovcagCg+55fI0n2KZmY/beyxbd/VmNMWw7CZmDjzvnvaHmJBNPnfK3oAzBUBEQlfp9AZOfu0hTAsD+N47DgDIiYALytAATnt55rUteeHbxQIH4sIKl5gxRxxBkUFPGur7h3/7RObz7Alt3PfKer74t+cCE10WisMXEV43YzTzjU7bXFipUKf/uTfunfPdG1319gc7j03GtDYwZaRjiuiKYBYDeP1eYxCBf7vLIAdx0NQO6pKSN+L1MDv2poBEhifvO46VW3bx0urO0EZ8zKwx9Kc1zwfm1aFQCPmJs8cxb9nmnOcT5NCHrFbxytrtOe+EObo+cXawYPE4tUhYb30qwVEzRnP/wg2Z7x4zx7Ywd3J7oUM5cd9xPPLqxswaIv1eAEhAXTznsd/s1tOvOeuymJy07zh29fZzz4J12f2LaCAHThlBe2Mq85y9QIedAfOFTDwTlLfEdiFBCY7J0Zuf42mVfSF1N/n48TP5z+eOi7T0s8dhe+QKl2IDt8Nd4XfPAucerNi8k3Wd3YHRnW/cbzx9aeXel7P3WFWdNlzkOnFghUvMmD4DsyPxRkX/9fNH8nJgffrEvXLWDjniv+/OG3EHTdbzOGiq4wPxQldH5AiX4Ec+e0J7TmO66Rln4l1vf5r6VPFRzwn7jOOO+Ws54ydOrL1/IqCfsW0NHDNrLL95aGlmmzmKdM6R4g17jeH2+WvyItn+/eJaPuqa+iA7qdPk5DnjSSWEG59aQbfbkR08rSNvv0P3GMnI5jruenFNzvZis57fdshk+tPKP57JmiCDHPrgzF3xQj5fdc2dj766kbe4cxOufO8heQ5mP6fPzQqX/SYFC4rj9h7Lkg07WLphRyY328yxLfz5Y0eFOqjfeuAk+tOayZflKRVBo/cfnX0wo1rq6e3XnAWtekOc7Ue7s+IvuO6pjK+mmHDxtM8/z3uNZ1/bkhEu/QERdiYjW+o5bu+x3PjUCvrTmhHmhTSRv3z8aMa1NfDE0k3s7OmLrLnUJRM5uc2i4K1V41Gsz58+poWDpnbwt6dWoKrcOd8RmE8GBCwcNKWDsW0NmX0gm37Kai7DEDNc0xwBn7BPNvzvMZ/54vT9J+akxgZ4wRdppBCc/xk4Y/+JpJVMGndzBBumUfz5Y0fx1gMn8fHjZvLvl9by4qptLFq3vajmAnDGAROBbIdRaHKYyWdOmpUzA//3H81fduf0uRNZsXkXT/lMhqZgASflvJ8xrQ2cOncCf3lyRWZp5S+/ad+8/VLJBCftO567X1qXk2uq2dWGCmkUM8e2cugeI/nDo8vp60+jqjy/cmtBx/mnT3RnhS92NIt/PJ0VSh0FTFYmo1sbePvBzqDjh+86KHCfU/abgAhc//hy1nd2M310M3d/7vhMnq9CzJ7QxpyJ7Vz7yDJeWLk1cx8Ck112NPHH846kP63caoTv9hZI5QLOO/jx4xyn952uEPfelWs+dHjBcr39EGey4pk/eyizlow/VDuIdx46hbXburnrxbWBSThNDprawY/PPpiu3jQ3P7Oq4ITjSvHlN83OfA6a8OznHYdOYcGaTp5YurlgXjVw2twpc8Zz78vrMhaEzJIGIZOa48IKl5gxtQbTob9ySxcLvn0aDalEXobdaaOaGdFUx7NffyPT3TxZ/3h6Ze7oPWSS16zxbew/eQTXP76cdFpz7NyF7O7gjOCvOOdgLjhuJm0NKd50xQOAkyKjGH6t48g9RxXYM/d6Xlr+z52yd2Ak25sOmMiIpjquvGdRwfNMGdnEqfuND/ztQ0dPZ+uuXn5+r3N8oY7pnCOm0dndx58MH5enuVzuTgQM4oLjZrJ8004+++dn2P8bdwLw75eCTX17jWtlzsR2rnrgVbp6+zPCC4rb3j1+8K4DefW/35QTrWQyuaOJ0+dO4Jf3vwqQM5M8DBHh48fPZNG67bz5Jw/y2T8/AxTulGZPaGffie38+N8LMwuqdfemA5+hx8Wnz2af8W1c89BSevrSdLsd34iAJYs9Tt53XEYr9XwthQSrySlzxjN9dDM/vnthxswVpiUcseco9pvUzvfvfCVTn1JyiJXCecfM4NFLTuJbZ+7Hf4e8Wx7/dcgUxrQ28KN/v0JdkTKdsf9EdvT056wPBFZzGVZ4KrW55rjnczl4WgffOnM/GuuSHLv3WG58aiX/et4ZzX3h1H0yq+yNaK7jns8fz3nH7MmNT6/MNb+EmMUAPnrMnixct53bXliTMQlBdlZyGCOa6/jiadnRVZScRYmE8MAXTwDgkGkdOR1nGN8+az+SCeGgAHMVOGvBnHfMnty9YB3vu/oxFgWkqfnrx4/OpAvxc/j0UZw4e1xG8ylkrjt0j5EcPXM0P793UeaZ9fSnOWDKiJxQYD8n7TuOA6eM4JbnVhfNgyYifPLEvXht0y5+ed+rOWUpZns3z1FMK/zsyVn/2aiWcI3F5Iz9J+YNPsLmYXzmpL3YuKOHq+5/lbtfWsvt89eEjqwBLnnTbJZs2MFV9y+m2/PPhQgkEeH3Hz0yZ5s5670QqaSTOPWl1du47tHlReuSTAifP3UfNmzv5ma3Y46rQxZxsil84KjpmawaYTTVJ7ng+Jk8vHhjpg/443lHBu579MzRzJ7Qxi/uW0xvfzrT50SxPlQaK1xiwjMPmFqDN4r48pv2zTjqLjplb7bu6uXCPz4FOLPsTUSEi0/fl0OmdXDpTfMzE9HCFjMCePMBk9hnfBtf+ttzXOWOYn/xvkNoaywuXMAZye893nnxvxJgSgpi6qhm/vrxo7jinIOL7+xy6B6jmP/NU0NnCX/QDcN9cNEGTv7h/TxvzEl5/Msn5c0t8PPVM7LlDxN6F52yNxu293DOrx515uz09Be1u4sIV3/wcN+2wvu/af+JnDh7HD+9Z2HOvI+omksU9h7fxhn7O2bKsCWn/SQTwu/PPYK2xuw9CgspP23uRE7edzy/uG8x517rmCmDknOaHL/PON58wET+798Luc911BcLWz9k2kh+86HDmDGmhT9+NLhTDeLMgybxhr3G8PhSx+xc7FkeN2sscya28707Xo60fzX54FF7MHdyOw8vdqwIB08NnrMiIvy/U/Zm8fodXPPQktDMHHFjhUtMeI5NM1wzaGLavhPb+cTxMzPfvXkqJsmE8OOzD6ahLsn7f/0YC9ZsY2dPf9E4/CvOOThnNH3a3ImRy59MCP/69DE88/VTOO/YGcUPcDls+qiCWkQhikWWtTXWcc2Hsx245wT/wqn7hGoVHjPGtvLjsw/i9LkTGBfS2R42fRTvOHQKz7y2hS/89TkeW7IpLzw8iLFtDdz3heMzDve7LzoudP9LTp9NWsmZc1PsHpTKj84+iCvOObhoBJqfA6Z08Pw3TuVrb57DnInteRmt/XzjrXNyRvhRfAjfffv+7DG6mesfd0yQUep+4uzx/Ofzx3P0XmOK7ushIvyPG8oL+QlB/SQSwg/edWDmeyUF/kBJJRP86N3OoG3/ySNC6/LGOeM5ed9xfP+OV7hhnrOMQ9SU/JWkdu7eMMN7MXf1ZE1SngPTH1HzhVP34ZZPvYH/e/eBBR2vU0c18/tzjyCtcNqPHuCBhRvy5hn42WdCGz9xtQhvxnoppJKJohFM1eKEfcbxyCUn5mwLihArxJkHTebn7zu0qEnpO2fNZcbYFm59znFUHxHBdwROmo+fv+9Qll5+RtEIolnj27j+vNcBjunl82/cO+NbqxR1yQRvPXBS2aPvc9+wJ//6zDFFO9gpI5u56gOHZb5ffPrskL0d2hrr+MO5Rxrfo5lQy2FyRxO3f/YYPnvyrMwS22HsO7Gd75w1l9aGVOCEyMFkr3GtPPGVk3PuXRAiwvfecSBjWuv5n9udBc3CfK1xIXGtnzyUOOyww3TevHnFdyyBQ799V2b50BljW7j2w0dw54tr+fYtL/Lvi45lr3HlvbirtuzivN/NY/6qbRwzawy/L/KigeMzEaFoxNBQIJ1WPnX909z6/Gqu/sBhnDwn2JE/EHr701z8t+f521MruPYjR3Dc3vEk9tvR3Ue/ak4CxaHKwrWdjGqpzwu1DWPD9m5Wbt7FgWUMfOKmrz89KHNDKsnKLbt4/eX/AeDmT74+Jx9apRCRJ1X1sMDfrHCJR7gc9K07M+Gvfh780gklm478LFrXSWtDXVF/w3BEVVm+aWdOUsA4WLVlV0kT5CyWWqOnL83i9dvZd2LhCbQDIUy4xKeP7uaErWldCft6uZrPcEBEYhcsgBUsliFPfSoRm2AphhUuMRGU9rw+mWDOpHZG1Ygfw2KxWOJi2AoXETkN+DGQBK5W1curef2+tHL+sTN456FTmDa6OXRymcVisQw3hrbHqgAikgR+BpwOzAHOEZE51SxD2k1zPWt8mxUsFotlt2O4ai5HAItU9VUAEfkTcCbwYiUvsq2rlxdXbaO7L013bz+rt3Zxy3OrmDm2ld7+8NUYLRaLZTgzXIXLZMBcCGUFEH1qb0ReXtPJ2Vc9mrf9iaVOxlJvhrvFYrHsbgxX4RI0cyzHwy4i5wPnA0ybNi1g9+LsM6GNP553JA2pBA2pJA2pBGl1TGKjW+ojzR63WCyW4chwFS4rgKnG9ynAKnMHVb0KuAqceS7lXKS9sS6zToXFYrFYsgxLhz7wBDBLRPYUkXrgbODmQS6TxWKx7DYMS81FVftE5JPAHTihyL9R1fmDXCyLxWLZbRiWwgVAVf8F/Guwy2GxWCy7I8PVLGaxWCyWQcQKF4vFYrFUHCtcLBaLxVJxrHCxWCwWS8Wx67kAIrIeWDaAU4wBNhTdq/YZLvUAW5daZLjUA2xdPPZQ1cDV9KxwqQAiMq/QgjlDieFSD7B1qUWGSz3A1iUK1ixmsVgslopjhYvFYrFYKo4VLpXhqsEuQIUYLvUAW5daZLjUA2xdimJ9LhaLxWKpOFZzsVgsFkvFscIlAiIStD7MkGQ41cViiZPh1FYGoy5WuESjbrALUEGGzTMXkTHu/yG9nrSIHCYi4wa7HJVAREYYn4d652zb/VC64FBCRI4Skb8A3xeROUO5ExORI0TkD8B3RWR/ERmSz14cmkXkeuAmAFXtH+RilYWI7CciDwOXAh2DXJwBISJHishNwNUi8hERadAh6tC17b4yDMkOphq4I8mf4qTt3wB8BviI+9uQGJG5HXFCRC4FrgZuw1lm4ULgwEEtXJmow0736xgRuQBgiArLzwB/V9W3qOorMHTeLRMROQD4GfBX4C/AicBeg1qoMhkO7R6c9jDY7X4oNshqcQCwUFWvAX4A3AicKSJ7q6oOhRfN7YjTwGvAh1T1OuAyYA+cRdSGHG6jmQisBc4FLhCRDlVNDyUBIyKjAcXpyBCRt4nIFKDJ/V7z75fBYcAiVf09cBfQCCz3fhxidTkQeGUot3sAt92vYBDb/ZBpjHEjIu8UkWtE5CJ30zPAoSIyU1V34CydPA/4GDgd9+CUtDhmXdzG8EfgWddUsRHoBCYObimjYdTlsyKSVNW0qq4GpgNLgfuAi93nlB7MsoZh1OP/uZt2AMcCJ7pmi48B3wF+BEPn/XI33QK8TUQuA54HpgBXiMiXYEjU5asicqq76WngsCHc7s26/BF4ZrDa/W4vXEQkJSL/C3wax4Z/gYh8VFU34Kj4n3J33QL8G2h2R841R1BdcEYuXarar6rdIlKH0/hfHsyyFiOgLhcCH3B/mw28qqorcEbKnwD+IiINbv1qhoB6fEJEzlfVLuBaHHPSHap6GvAVYK6InD54JS5MgbZyvqquAw7GMb18WVVfB/wWeIOIHDVoBQ7BrcvlOHV5BfimiHwAp53/AcccBkOn3fvr8l5wNJjBavfDdpnjqKhqn4iMxWkUD7imiXr35z8D/yMiJ6vqv0VkIzAZ2DpY5Q2jQF2afbvNAdaq6isi0gbsq6qPV72wRShQlxb35+XALBH5J7A3jvbSp6rdg1TcghSoR4P780+AzwLj3H1XisiDQE1qYGFtRVVfdIX+He7uTwLrgJp7JpCpy1TgIlV9QkS24miNa4B/At8ZYu0+qC7rcAZfAPtS5Xa/W2ouIvIBETlORNpFpAHYBnxcRL4G/BA4RkS+gJOG/1rgxyKyF3ASIGSFz6AToS6vF5GLRWSUe8hoYKeIfAh4GNi/VuzIEepytIh8EWcEthRYAhyqqm8BporIoYNVdpMI9ThKRL7iml0+BnxQRA4SJzjhZJy61QQRn8lXRaQRuBn4hvs+nQ3sB2wctML78NWlHdiFayZS1Tvc70cDq4A/AT8aIu0+qC47gWMNbavq7X63Sf/i3sgJOHbINLAYZyT8CZyX5r3AO4D34TgkvwX8Q1X/6HZoewOzgfNU9aXq1yBLmXX5u6peLyLfBb6EY7b4kao+V/UKGJRRl2/iOFlvV9WtxnlGmN+rTZnP5CZVvU5EzsYJINkPRyuYX/0aZCnzmfwdR7j8EhiLYxX5tKq+WO3ym4TU5XwcAXg68CrOyH6p+/8CVV0oIhfjRL3VersPqssSHCvFx1V1kYj8D/AFqtnuVXXY/wFJ9//ewB/czykce/eN7ve5wC+NY74O/Nr4Xj/Y9RhgXa52P78BePdg12OAdbnK/ZwAEkO4Hlcb32Ww6zGAulwK/Mr9XA9MGOx6FKnLlcDv3e8H4gy2/sv9fg2OcPHOUevtvlhdPu5+Prra7X5Y+1xEJIUzQkyKyL+AdqAfMnbKTwKrReT1OOp+v4icqKr/wVmd7RXvXKraU/UKGFSgLgvdfR8clAoYVKAui919B9U3Ualn4u4/qCaEAdZlNLDI3bcHx28xaESsyyoROUFV7xGR59UJZW/EMX/d7p1riLT7sLrc4e77cLXLPmx9LiJyHI5TcSTOi/9toBc4QUSOgEyD/hbwVVV9HsfWepGIPA2MAn4+GGX3Y+tSe3UZLvWA3bIuaeAbOBoXQEpEzgNewAkRX1ULfsghX5fBVvdiVCOPAd5vfL8SNzQXeNLdlsCxYf4NV5UHZgGzB7v8ti61XZfhUo/dvC434KwBD3AWcPRgl3841WXYai44Ev8GyeYFegiYpqq/xVExP6WO1J8C9KrqGgBVXaiqCwalxIWxdam9ugyXesDuW5d+VV0GoKr/0EEwHRVhSNdl2AoXVd2pqt2aTWp4CrDe/fxhYF8RuQW4Huch1iy2LrXHcKkH7NZ1eQpqNz3NUK/LsHboQyYduwLjcUIlwUmD8GWcqJclqrpykIpXErYutcdwqQfsvnVR15ZUqwzVugxbzcUgjbMuwwbgAFfSfw1Iq+qDQ6WxuNi61B7DpR5g61KrDMm67BaTKEXkdTizUh8GrlHVXw9ykcrG1qX2GC71AFuXWmUo1mV3ES5TgPcDP9QazD9VCrYutcdwqQfYutQqQ7Euu4VwsVgsFkt12R18LhaLxWKpMla4WCwWi6XiWOFisVgslopjhYvFYrFYKo4VLhaLxWKpOFa4WCyDgIh0iMgn3M+TROSvg10mi6WS2FBki2UQEJHpwC2qOnewy2KxxMGwzy1msdQolwMzReQZnEXD9lXVueKscX4WkMTJG/UDnNUd3w90A29S1U0iMhNndcixOOuln1eDGYotuzHWLGaxDA4XA4tV9SCctc1N5gLvAY4ALgN2qurBwCPAB9x9rgI+paqHAp/HWevDYqkZrOZisdQe96hqJ9ApIluBf7rbn8dJXNiKsyb6X4wM6w3VL6bFUhgrXCyW2sPMHZU2vqdx2mwC2OJqPRZLTWLNYhbL4NAJtJVzoKpuA5aIyDvBWSBKRA6sZOEsloFihYvFMgio6kbgIRF5AfheGad4L3CuiDwLzAfOrGT5LJaBYkORLRaLxVJxrOZisVgslopjhYvFYrFYKo4VLhaLxWKpOFa4WCwWi6XiWOFisVgslopjhYvFYrFYKo4VLhaLxWKpOFa4WCwWi6Xi/H8Ii5ScUAVgIAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Start a new model instance\n", + "model = GR4JCN()\n", + "\n", + "# Accessing the future data we just created -- notice it is now the \"future_dataset.nc\" instead of the \"reference_dataset.nc\"\n", + "# Model configuration parameters\n", + "model(ts='future_dataset.nc',\n", + " start_date=future_start_day+dt.timedelta(days=1), # Now running the future period!\n", + " end_date=future_end_day-dt.timedelta(days=1),\n", + " **common_run_parameters, \n", + ")\n", + "\n", + "# And plot the results! We can compare both figures and see the changes between both.\n", + "model.hydrograph.q_sim.plot()\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "cc451f43-c827-439d-adef-713c02965575", + "metadata": {}, + "source": [ + "## End\n", + "\n", + "We could now download the data to a netcdf file or other medium, and we could analyze the data either locally or on PAVICS-Hydro by appending other blocks of code below where the analysis is performed. " + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.8.13" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/docs/source/notebooks/paper/Qobs_574_HYSETS.nc b/docs/source/notebooks/paper/Qobs_574_HYSETS.nc new file mode 100644 index 00000000..9e5b12e5 Binary files /dev/null and b/docs/source/notebooks/paper/Qobs_574_HYSETS.nc differ diff --git a/docs/source/notebooks/paper/shapefile_basin_574_HYSETS.zip b/docs/source/notebooks/paper/shapefile_basin_574_HYSETS.zip new file mode 100644 index 00000000..ead0a008 Binary files /dev/null and b/docs/source/notebooks/paper/shapefile_basin_574_HYSETS.zip differ