diff --git a/.gitignore b/.gitignore index 43871650ed7c45caf2bef3a7073e5a9fe788cbff..1ee30200347a98067da3c2c40d136c735823d755 100644 --- a/.gitignore +++ b/.gitignore @@ -2,3 +2,5 @@ *.nc *.grb *.grib +notebooks/dask-worker-space/* +*png diff --git a/notebooks/use-case_climate_station_data.ipynb b/notebooks/use-case_climate_station_data.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..1bd6684f9f67faefd26b4b2b3de328ffd7fa208b --- /dev/null +++ b/notebooks/use-case_climate_station_data.ipynb @@ -0,0 +1,400 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEKCAYAAAAVaT4rAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABQNElEQVR4nO29e5gcZZn+fz8zmcxMkkkmmRzIiSQk4RCIBJKAIAGWFcQjymEVdZFLFP3iqsiKoq4K6nqxKheyHgFBYL/C76t4WhBWQYGwyilAIOFsIIRJQk6TyWEyM8lk3t8fTz/U29V17qru6u7nc11z9UxNd3V1ddVdd93v+z4vGWOgKIqiNA5N1d4ARVEUpbKo8CuKojQYKvyKoigNhgq/oihKg6HCryiK0mCo8CuKojQYmQk/Ed1IRJuJaLW17HIiWk9EKws/78jq/RVFURRvsnT8NwE43WP51caYRYWfuzJ8f0VRFMWDzITfGLMcQE9W61cURVGSMaIK7/kvRHQegBUA/tUYs93rSUR0IYALAWD06NGLDz300ApuoqIoSg0zMAA88wweB7YaYya5/01ZlmwgotkA7jTGHFH4ewqArQAMgG8CmGqM+WjYepYsWWJWrFiR2XYqiqLUFU88ASxeDAIeN8Yscf+7or16jDGbjDH7jTHDAK4HcEwl319RFKUh2LMn8N8VFX4immr9+T4Aq/2eqyiKoiQkRPgzy/iJ6DYAJwOYSETdAL4O4GQiWgSOetYC+ERW768oitKwVEv4jTHneiy+Ia3179u3D93d3RgYGEhrlbmmra0NM2bMQEtLS7U3RVGUarFnD9DTA8yYEf68AKrRqycVuru70dHRgdmzZ4OIqr05mWKMwbZt29Dd3Y05c+ZUe3MURakWV10F/OhHwOuvBz8vTxl/mgwMDKCrq6vuRR8AiAhdXV0Nc3ejKIoPmzbxT39/8PPqVfgBNIToC430WRVF8UHM37Ztwc+rZ+FXFEVpKET4t24Nfp4Kfzb09vbixz/+cazX/PSnP8Utt9yS0RYpilL3xHH8o0b5/luFPyFJhP+Tn/wkzjvvvIy2SFGUuieO4w8Q/prt1VNtLrvsMqxZswaLFi1CS0sLRo0ahSlTpmDlypU488wzsXDhQlxzzTXo7+/H7373O8ydOxeXX345xowZg89//vM4+eSTceyxx+K+++5Db28vbrjhBixbtqzaH0tRlDyTkuOvC+G/+H8uxsrXV6a6zkUHLML3T/++7/+vvPJKrF69GitXrsT999+P9773vXjuuecwYcIEHHTQQfjYxz6GRx99FNdccw1+8IMf4PvfL13X0NAQHn30Udx111244oorcO+996b6GRRFqTNScvwa9aTE0qVLMXXqVLS2tmLu3Lk47bTTAAALFy7E2rVrPV9z5plnAgAWL17s+xxFUZQ3UMfvEOTMK0Vra+sbvzc1Nb3xd1NTE4aGhgJf09zc7PscRVGUN1DHX106Ojqwa9euam+GoiiNhDr+6tLV1YW3vOUtOOKII9De3o4pU6ZUe5MURal34jj+qVN9/53pRCxp4TURy3PPPYfDDjusSltUHRrxMyuKYjFxIrv92bOBV17xf968ecCxx4JuvbX6E7EoiqIoZaAZv6IoSoMxMAA0NQG7dwODg/7PU+FXFEWpA4aGgP37new+qIFXhV9RFKUOkJhHJmHxE/59+/hHhV+pKps3A7fdVu2tUJTaRoR/+nR+9Mv5pVa/Cr9SVX7xC+CDHwR6e6u9JYpSu0R1/FKSWYU/fZJU52xY+vqKHxVFiU9Ux6/Cnx0q/DGQAzZsujhFqTbGAFdcAaxbV+0tKcUt/Or4K49dlnnp0qU46aST8E//9E84+OCDcdlll+EXv/gFjjnmGCxcuBBr1qwBANxxxx049thjcdRRR+Gtb30rNm3aBAD4zGc+g2984xsAgD/+8Y848cQTMTw8XLXPljoq/EqtsGEDcPnlHE/mDTmPxo0Dxowpy/HXR8mGiy8GVq5Md52LFgEepZSFJGWZTzjhBDz88MMgIvzsZz/Dd77zHVx11VW48sorsXTpUixbtgyf+cxncNddd6GpqY6uySr8Sq0gcWTQqNhqIedRW5szgteLhhH+HCBlmQGUlGW+7777AADd3d14//vfj40bN2Lv3r2YM2cOAGDUqFG4/vrrceKJJ+Lqq6/G3Llzq/MhskKFX6kVakX4u7rU8Qc580oRpSzzpz/9aVxyySV4z3veg/vvvx+XX375G69ZtWoVurq6sGHDhopud0VQ4VdqhVoR/jIdfx3lCZUlSVnmHTt2YHqhYebmm29+Y/mrr76Kq666Ck8++STuvvtuPPLII6lua9URwVfhV/KOCP+6dTxKNk+k6PhV+BNil2W+9NJLI73m8ssvxznnnINly5Zh4sSJAABjDC644AJ873vfw7Rp03DDDTfgYx/7GAbkS64H1PErtYII/759wPr11d0WN5rx54Nbb73Vc/n999//xu8nn3wyTj75ZADAGWecgTPOOKPk+fZcu4sXL8aqVatS3c6qo8Kv1Ar2WJNXXgEOPLB62+LG7fh37OALVEtL8fPU8Su5QIVfqRVENIH85fxuxw8APT2lz1PhV3KBHLD2SaUoecTt+POE2/ED3jn/nj3AiBGldwIWNS38tTB7WFrU9GdVx6/UCiL8U6fmW/jF8Xvl/CElmYEaFv62tjZs27attgUxIsYYbNu2DW1tbdXelGSo8Cu1Ql8fMHIkMH8+8PLL5a3r3/+dBfptbwO+/nXg2WfLW9/AANDczG4+zPGHCH/NNu7OmDED3d3d2LJlS7U3pSK0tbVhhlTlqzVU+JVaoa+PRXPOHOCee4Kfu2cPZ+x+5+UDD/Djpk3At74FfOMbwLvfDXz5y8Cb3xx/2wYG2O0DZTv+mhX+lpaWN0a+KjlHhV+pFfbsAUaPZuHfsKFYbN38x38A114LvP669/+7u4GTTgJ+/Wt25j/6EfCDHwAnngisWQPMnBlv2+xtKdPx12zUo9QQKvxKrdDX5wg/ALz6qv9zX3yR3bzXmBtjgNdec8R94kSOe+64g7tgPvlk/G2zhb+9ncU9b8JPRDcS0WYiWu3xv88TkSGiiVm9v5ITjFHhV2oHEf6DDuK/gxp4pbyKV9yyYwdPiO529QsW8ONzz8XfNvfdx4EHAn//e+nzquz4bwJwunshEc0EcCqAHBa8VlJn3z5ASkyr8CtJWbECOP98YHAw2/dxO/4owu/Vl767mx/d+f+4cVxPP0lDr1v4Fy8GHn+89HnVFH5jzHIAHnsEVwP4AoD6746jFN8Gq/ArSbn7buDmm4Gf/CTb9xHhnzoVaG31F35jgh3/a6/xo1eOf9hh6Qn/+vWlbQx5y/iJ6D0A1htjnorw3AuJaAURrWiUnjt1iQq/kgbbt/PjN7+Z7dzN0qunqQmYNctf+HfudAYkxhX+BQs46vGabKm/3/+i4Bb+JUv40e368yT8RDQKwFcAfC3K840x1xljlhhjlkyaNCnbjVOyQ4VfSYPeXha97duBK6/M7n2kVw/AcY+f8Nvl072Ev7ubLx6FOTqKWLCALzBycbC56Sbg6KO956d2C/9RRwFEHIO5P0NehB/AXABzADxFRGsBzADwBBEdUMFtUCqNCr+SBtu3A/PmAR/+MHDNNd6imQYS9QDlCf9rr7Hoj/DoMS8NvF7OfutWbsfYvLn0fwMDHD8JY8YAhx6ab8dvjFlljJlsjJltjJkNoBvA0cYYn06wSl0gwj92rAq/kpzeXqCzk6MeY7gPfRbYwj97Njfc7t5d+jxb+L0ad+2unG6CevbI+eJ1MfEaU7BkSbHjN6bq3TlvA/AQgEOIqJuILsjqvZQcIwfy+PFapO3VV7kftxKf3l4+hmbNAo48kgdApY2Ipgi/DJLyEnYR/vHj/aMevxG9XV3A5Mnejl96LXn1z/cS/sWLgY0bne0ZHOTPUcVePecaY6YaY1qMMTOMMTe4/j/bGOMzhYxSN9jC3+iO/8c/Bt7//mpvRW2yfTs7fgDo6ABizn4Xif7+YtGU9/NqTN6wge9iZ84sFX734C0vFizwFv4kjh9w4p4IJZkBHbmrZI1b+BugqJ4vfX28D/I2pV8tII4f4GzbK34pF2lQFccv7yc9imw2bACmTWP37hbp7dtZgINqa4nwu88HOV+iOv5Fi7gRWeIeFX4lF9jCbwywd291t6ea6AjmZOzfzyNhxYF7Cf/jjwNXX13e+4hoivCHOX4RfncUJIO3ghz/YYfxZ9q4sXi5RD1RHf/o0bwudfxKrhCxmzCBHxtZ9HRCmmTs3MmPQY7/lluAL36xvPfxc/xewr9+vb/jD+rDL/j17Inr+AGngVfaKAAVfqXK2I4faGzhFzenwh8PiVqCHP+uXVweZN++5O/jFn55P3fUI6N2p01jQ9PTUxzZ+JVrsPHr2eOX8Q8N8Z2Pl/AvXszF4tavV+FXcoIIvQq/Ov6kiOO2HX9fX/HIV2nsLWffuoV/3Lji9xe2beMLjDj+oSHnrgRgx9/c7D14S5gyhT+P2/H79eqxZ99yI7X9//pXFX4lJ6jjd1DhT4bb8Xd08KM9ulWENw3hF9FsbuaeO27HL10nRfiB4pz/tdf4f83N/u9FxNn8888XL/dz/EHCf/TRfJG6914VfiUnaMbvoMKfDC/HDxTHPeL4vUodRMXduCvv6Xb8XsJvC3VQH36bCROK7xQA/4w/SPibm4FTTuEZw9wXLx9U+JVskQNW3FoWwn/ttcB3v5v+etNGhT8ZXhk/kL7wu6MeeU8/xz99umNobOEP68MvtLWVng92rx673SBI+AHgrW/lAYKrVvHfKvxKVRkY4Mmr5WTKQvh/9SvghhvCn1dtVPiTIY47ivC79600tEbBS/iDHP/UqaWOP8rgLaG9vXT2Lvl7YKD4s0QRfgD47//mRxV+papIF7T2dv47C+EfHCyunZJXGln49+8Hbr/duxRxGNu38yAlyfa9hF8iE9vxP/YYC/ATT0R7nziOv6uLC6a5hb+nh7/nKFFPe7u34ycqXicQLvzz5/NnfeYZ/luFX6kqlRD+gQF2fFmM5kyTWhH+lSuBY45JtyzCXXcB55wDPPpo/NdKgTYRRBF+2T5jvB2/9Kd3ly32Q4RfjlXA3/FPm+b8H3Aad6P04Re8op6BAa7jAxTn/GHCT+S4fvdn8ECFX8mWSjl+oHQUZNqsWQOcdZb35NpRqBXhf/BBdssvvpjeOqXb4o4d8V9rl2sAHOcvF/rBQe5SCRQ7fvl/1Nmu+vr4OG2yZNHP8YvwjxjBPWrEnUvxuFmzwt/PL+qZPp1/j+P4AeDUU53nNAVLuwp/pdmzh/vwNkqVRrfwZyF6lRL+5cuB3/wGePnlZK+vFeHftIkfvUaPJkW6LSb57HaBNqA06rHvTMoRfrsypzB+PK/THhhmCz9QPHr3qadYdA8/PPz92tr42LXjr8FBR/jjOH6Ae/YAoTEPoMJfeTZu5EkWXnih2ltSGSoV9QDJc/7LLnMaxYKQbU8q3LUi/DIJSJpTnpYj/G7HHyT89vrl/151772wa/EL7no9+/fzHLd+wr9yJU+OEhK1AHCeI8eFMeU5/ilTgDe9SYU/l0gWmDQuqDXyHvX09gLf+Q7w61+HP1e2Pcln2L/fiSPyLvzi+NMSfmOydfx2X3gvx9/dXdpf3gsv4XfX69myhb9LW/ilbAPAjv/II8PfCygVfilgKCN+4zp+ALjkEuCf/zn0rVX4K43khY0i/P39fLC2tHAemqXwJ3H8f/sbC1OU7Lkcxy/bmPT1lSTtqGfzZkc403D8I0fy8SROP8zxA9Fcf5Djl/N2/Xp+9HL8PT3AunVcKjkKIuJyXMkxMnp06QQvUYX/Ix8Bvv3t0LdW4a80cgA1yghWu6KgV/e1NCjH8T/4ID9GcYTlOH77Ql8rwp+W47fLEqTh+AFu4I2S8UtPoKjC745J3I7fq9eOCP/TT/PfcR2/HE+2uHd1JXP8EVHhrzSNGvUA2Qi/5KJAMsf/v//Lj1GEv5x6+rUi/Mak7/jLEf6BAb6w244fKK7Q6ef4+/q4d01ra7QG3iiOf906fnQL/44dTrfRqI7fHfXY4j5xYjLHHxGPKeCVTGm0qCdr4R8acoa2x3X8AwNOv/Kso55aEX6ZJQxI1/GPGsXuO+5nd5drELyEf+TIUsff2cl3B1GE369XD1Ds+FtbgUmTnOdI2Yb77+cG1ilTwt8L8I96ZGCYbWTU8dc44vgbKeoRZ5OF8NvZeVzH/9hj3KA2bZpGPYK4fSBdx3/IISyqcT+7u0CbYAu/fHdTppQK/5gxXPs+rYxfyjFIhAQ4o3cfeCC62weCox4vx9/UxO1kKaDCX2nU8ae7fhF+Ee84Rbok33/729XxCyL8c+ak6/gPPZRdf5aOf8qU0sZdEf5XXgl/by/hb2/nOwm5AK1bVzoqV4R/9+7o+b6sG/COerwy/ra24gtOGajwVxrN+IOf/+ij8QpryX6cPZsf48Q9//u/LApz5vAFxL578CINx9/RURvCf8QRfKyWOzF8fz9XjUwq/O4CbcKYMcW9ekaN4n3r5fgPO4zjwLCxM17CT1Q8eve114ADDyx+jkQ9QDzHHxT1TJzI+8q+G0gp5gFU+CuPOv7g5595JvD1r0dfv5wsc+bwY1Th37+fZyxatsyZaSmsNk0ajn/ChNoQ/sMPZ7F0TyQel5de4vWI8Mctmyznizvqcffq6egovbDYjh8ozfkfeYSPm23b+HgYHPQe/CT1eoaGOE70c/xAelGPu/ibCn+N00gZvzF8MsUR/i1bnHonUXALf9Scf9UqjoZOOIFnWQLC4540evWMH59v4ZdRuyKW5cY90qMnC8fvFv7Ro70d//z5PFmJO+dfvRpYuxZ48knvSVgEcfwbNnB5BbfjF5Fua+P3ikpYrx5Ahb9uaCTHL6IcVfgHBrix9dVXo7+H7Me4jl+6cS5b5gh/WANvGlFPLTj+8eOdAUrlNvA+/zzHJfPnJ2vcjZLx79zJ36H7jkKEf+RIYN68Uscv38nzz3uXZBbE8ftV3hw7li8sCxfGa3wN69UDOPtfhb/GaSThd3dBa28PPvHFcXd3O+UNwpCT5YAD+AQXx3/nncDixcXFtWxeeIHFZNYsJ+oJc/xpjNwV4bdnV8oTmzZxI6k4zjQc/+zZ/N0ndfzt7SyGNvaE67bjl/Xv3cs/Ut4haH7bMOEXx+8n/ERsPI4/Pt5nq6Lj1378lWTvXucAq0TUc9llLGY/+Un27+WFW/hHjQr+3HJbL1mq+5baCxHU9naucSKO/9preQKOzZudolc2mzc7/a0r7fj37+cL0siR8deTNSL80k89Dcd/6KH8e1Lhd+f7gCPofX0s/NOnF0c98ijPmzCh9MIe1/F7Dd4S/vpX572i4s747TtkOSbV8dcBdl3vSjj+++/nBqxqIQd01KjHPjGjxj327fG0aXzB6OsD7r2Xl/s1Tm7e7Ex4IY4/qvCX27ibdB2VQIRfooZyHf+6dU6Pq6TdOd0xD1Bck99u3B0c5AurxEAixkHTHD73XLjwi+MfN84RZZvJkyNVxSyipYXvFtyNu62tpXP5qvDXMCJCkydXRvg3bQrvopglXlFPFMcPRBd++2QRx3/PPc5yexCMjS38URt303L8QP6Fv7WV90s5jt8YvpiKcGfh+G3hF9Hes6dU+Nva/IV/wwbnTtFLvDs7+WLy3HPRZtaKClHxBck+X0aOLN7/Kvw1jDj+qVMrE/Vs3uyUeq0GXsIvMYcX5Tj+tjbH8f/+987/ozj+uFFPvTr+gQHeBxKBTZxYnuMfHOTvWvavCH+c9g0/x28Lv924C7B793L8/f3F721fCJ58kh/9HD/ARdiixI9xsKdfdHeGsAdxqfDXMCJC06axICeZeDoqfX18kuXB8dslGwD/i54If1NTsqhn6lRex+9/z711AG/hHxriOwERfnFYQY5fuqYGbX8Qsi9ExPIo/NKVU/bLpEnlCb9cSG3hD7rwexHm+Ht7+fuI4viNKX5v+3t8/HF+9GvcBXhfpOn4geK7YDlGpO1n7lxn8nQV/hpGHL90lctSlOUkzoPw244f8BdOiXrmz0+e8QO8n88/n3/3Ev5t21gE7EJbY8cGO37bHSYV/ra2YnFKg/5+7o+eBjJ4y3b85UQ9XsIPxPvsYY7/9df50RZ+L8cvx6D9PQ4MsIMfMSJY+O0LTxbCb0c9ra1OWYZly/guY/t2Ff6axnb8QLZxjwh/3qIeINjxy3ylSTN+gBvNzj6bH72EX1ysOFuAG+2CHL+9zUmjnra2ZOIXxI03AkcfHa3IXBhu4c/C8QPRP/vwMH8nQY5fsnlp3JX1e0U9QPH3ODDAr5s3zzlfghw/kH3UY3dbPekkNih//asKf00jjv+AA/gxywbeWnX8Y8dyL5B166Jlwe6MHwD+4R94PfaUeDbuSAMId/yyzUkLzQUJfzmuesMGji9kZqhy8BL+rVuTjzkoV/h37uT3DurVI+M2kjr+tjanuylQHcfvV4/nmGM49nnggdoRfiK6kYg2E9Fqa9k3iehpIlpJRH8iomlB66g7tm9nZykHY5bCLydxLQn/jh3OoKr+/mhu0456Zs1iAfjQh3hZFsKfdOSt3Ma7xe+JJ3g7ok4I7kbMhEQe5SD7xY56Bgbi19cRyhV+t3jbyDJb+JM4/rY2HtwF8N2m19gK+8KTddRji3t7O4v/8uWVE34iaiOiSR7LJxNRlC24CcDprmXfNca8yRizCMCdAL4WY1trn54edg9e7iNt5CQeHi6/wmJSkgj/uHEs4EC0uMduEBs7li94Mtl0V1d04Y8a9UyYwA477j71c/wvvMCuNsnsYYDTLpKG8G/axAIq31O5g7jKFX535wAbd9Qzdqy345dlURz/6NHeZY9lnAcAzJgRbdujEhT1AMCJJ3L7w/79FXP8/wlgmcfyUwFcHbZiY8xyAD2uZbalGg0gp+PWM2L7dhYOd42OLBBxA6rn+pNEPeL4AS6gFcbgIIu+nLDt7c7vQY6/ubn4Fj6O4w/6DEHb6SX8cleT9FgQx59kvmE30odfKLdsQ7nCb8drbmTCdTvjdwt/W5tTOyeq8HvR3MyfQcY3pElQ1ANwzi8mo0LCf4Ix5jfuhcaYXwA4MekbEtG/E9FrAD6EAMdPRBcS0QoiWrElrQkhqo04fneNjiywhb9aDbyVcPx29U83QcI/aRLf2gthjr/cfvh+jl+O7aTHQtqO374Lqrbjd4/8djNmTHDjrh0RBUU9hxzCf/sJP8Dnbdoxj2yXX9QDAMcdxxceoGLCHzTVS+K2AWPMV4wxMwH8AsC/BDzvOmPMEmPMkkmTShKn2sTt+Csl/JV0/L29zmhZ94kbduL39rIAyzypUYXfz4UFCb8tcIDj+P0aMuWzSCmDuA5dTmoRoLwKf9qOv6XF+X7sAVZRCIp6ABZ2ubh4OX5b+IMc/7hx3CMsSPgPPhg46qho2x2HsKino4N7bclzUyJIwDcT0THuhUS0FEAaFvxWAGelsJ7awZ3x12PU84lPAGcVvlb3gJSojbtE3LMnasYfJPx9faWfXxy/zdixPLDLT4DdUU9Sx9/czNublvCnGfXYhesAZx+VI/xjxzrRW5pRD+D07JHf5bwKcvxewg+wuEpvOy/uuAP40Y+ibXccwqIegHN+oGLCfymAXxLR5UT07sLPFQB+WfhfbIjInqXgPQCeT7KemsSYYMd//vnAxz+e3vtt2uScaFGinrTKBHd3AytX8vrc84QGCf/wMAuFNKTNmpWO4wdKXf+WLaWOP6w0c1qOHyiuWWMPyU9CWo5/xw6+U7MrmY4dy449adSzY0dxQbO4g9eiRD2Ac1fR1OTU5Pdz/Pb31t/vLL/5ZuD//l//bWlt5fdJm7CoB+CcH6iM8BtjHgVwDDjyOb/wQwCONcaElnwkotsAPATgECLqJqILAFxJRKuJ6GkApwH4bNmfoFbo6+PeIH4Z/8MPAz/7GVfULJf9+/lklR4IYY5fJt+QyUnKYdcuPuG3bvXungZ4i+bu3Sz+SYQ/KOMHSoXfL+oB/Bt403L8QLHwl9O429/vfLflCv/y5XyxPuEEZxlRefV6xPELafbqARxh7+govquIG/UAfEGvRqQcFvUAwGmnAV/6EvCP/5ja2wbW4zfGbAYQYwLUotee67H4hiTrqgvkltyvV490P/vsZ7n7VpyZfNz09DhTxL34Yrjwv/oqi/Xddxef+EmQeWtfeolPLPukDRJ+cdrSZ3rWLF4mDb5+xHX8UogsqeNP2qsnTPiTOH5x+1KVdO/e4Br/y5YB73wnz9Pg5r77eD+++c3Fy8sp2+AWfnf7RhhhUY8t/IJMxrJ7d/F37D727DvSaiKFCyVm9Nqe1lbg299O9W115G6lsCeN9nIfu3bx0PGnnwauv76895J8X3ohhEU9ctF57LF473PHHaUjRkX4X3yx9EBuaeHb8SDhtx0/EO76wzJ+oFj4vco1AOGO392rJw3h37+/uN56XOSYkgFIdruOF08/DaxY4f2/v/wFeMtbSoVnyhRn5qm4uIW/qYnXn3bU4xb+KI5/aIjNkd9FpVLYJrCCFyIV/koh4jNhQmnUYwwfqB/4AJcb+Ld/8y8nHAURgKhRj4j1Y49Fz/r37wfe977S2b3cjt8+kKX+uJdouifVFmH2q6cvBEU9ksfb+9Jr8BZQWpP/wQc59xXSiHrs3i179vB2yf4ux/GL8AfFPcawIMqIbptt24CnnuJjz83JJ/Po4iQlIdzCD8SryR8W9Yjg28Iv6w8TfndX42ph34kE3b2mTCzhJ6ImIvKYfkYJxXb88uXa9d2Hh/kA/ta3WBBkBqkkuB1/mPCL4+/tBf7+92jvsXs3i789q5jMcwp4Cz/gL/xuxy+vC9v2uFGPn/C7Z+G68krOVYX+fmcgj/wdBy/Hb2fn5Th+GYAU1LNn717+vryE/4EH+NFL+M8+mx9//ev421eu8EeNetwNyF6O3323mTfhHxjIl+MnoluJaCwRjQbwLIAXiChRr56Gxnb8I0bwjxx8IrwdHY5LF+ecBDm540Y9QPS4RwTSjkbs9cQVfnGvbuEPE8SgqKejg8U6juOXz/PMM8V5f3+/M2E4EM/xSy3/tIU/juOXvvNewn/ffbxNS5eW/u+QQ4CFC4Hbb4+/fWk4fiL/3jReUc+oUXzu7NlTLPxExbNw5UX47agn6O41ZaI4/gWFUgvvBXAXgAMB/HOWG1WX2I4fKO7GJSI/ZkzxIJSkbN7M7kb6JUd1/CNGRBd+2WZb+GVZRwcLv91dTghz/BL1RBX+IMdPVDqIS4Tfqx+/bMfu3dy2sGePc9EU4Q8bi+CFrCNI+JP06hHhl5GnQY5fjqedO0v36V/+wg2/fg3DZ5/NPb7ijBXYu5ffp1zHb5fgcOOX8ct+dQ/I8pvmsJrI8bR7N7c75CjqaSGiFrDw/94Ysw+NVmMnDXp6WFi95gC1BVP+X67wT5rkHFRRM/4lS4BHH432Hl6OX9Zz1FG8/WvXlp5YY8Z495zxi3qiCH/QyesW/i1b+Pnuio8tLby/du4Enn22dLvcwh/H8btFxi38XV3lRT2TJ/M6ojh+oLgReNMm/rxeMY9w9tl81/Kbkgou/sixkIbw++En/PL53N+x3XUyb8IvF/EcOf5rAawFF1VbTkSzAKQw60ODsX07u31xL/ZBaEc9I0dyPFGu8NsFpaJEPaNHc1e+J59k5xFGkOOXIebd3aUH8qGHOtPJ2fT28meX56fh+AFvxz95sreLHDuWhT5I+EeM4ItEHIceJvwzZiSPekaN4v12wAHRhd+Oe2TcSJDwL1jAP3HiHnedHiFu1BMkhH5RjxzvXsKfN8cv75834TfG/KcxZrox5h2GeRVAwFGieCLlGgS/qIeIRdjOy+Mi4ia37lGinjFjOOPt7/cWZjdBwr94sbPMfSAfeSR3D3T3WpJyDe7XlZPxAyz8ds8gr8Fbwrhx/Hnszy8npC1CUSZjuegibqiX1wLewt/ZyeKY1PHLPpO+/H74Cf8DDxTXg/Hj7LN5kJdXG4EXaQh/mOP36tVjxztu4fea37bawi+fT+7e8hL1ENEUIrqBiO4u/L0AwEcy37JaYMsW4Mwzw7scAt7C5hX1AE7PhKSIuMlBFCXqGTOGJ30AosU9cmLbsY18jgULnIuOl/AD3KfcRgq0CVk7fi+kUNszzzh3BCL8tghFEa+//AX485+LP4Mt/P39vC0TJxYfC3GwJyJP6vhffJG/r7ABg2edxT3P7r472rZVQvi9evVI47v9f8Hez2FjBCpFjqOemwD8EYDMlvUigIsz2p7aYvly4Le/jSaUXv2KvaIeILrwr1njXWbBLfxRop4xY4C5c1lI3A28+/YB111XHAF5OX75HJ2dwEEH8e/uA3nRIn5cubJ4ufvCKNuedsYf5vh37GDhP+IIZ7uAYhGK4vj7+krr8LirlK5bx20xSYXfdvwi/H7jMPyE/9VXncFyQUiX0e7uaNtWragnzPHnPerJi+MHMNEY80sAwwBgjBkCUKUpnXLGunX8GGWwVV9fqfB7RT3yGEX4v/EN4MMfLl7W38/rixv1SL2TpUtLhf9Pf+Kqm/ZFRk7swUFn/fady/xCPT63Y5syhX+eeqp4ubs0AxGfBEGCODTEfdODTpauLt6ufftYFMMc//r1/L2+5S28zMvxt7eHi1dc4U/aq0cc/9Sp/D5+JSe8GneHh3kbZs8Of6+RI/l7jXJ3C1TG8c+dywZD7iKBYOHPc+OuRD05cvx9RNSFQk8eInozgIAZKxqIOMIvDaiCV8Yf1/Fv3lx6otv91KX/c9SMH+Csd/Xq4qkFpWyC/TntcQbyu30BE+H3OpCPPLJU+N1Rj7w2SPjt+Xb9kEFc27fz9g0OBgu/fNbjj3e2CyiNeqI6fqkJI59HXg9wJp+m4wf84x45njo6HMf/+ut8NxjF8QN8EY1atydI+KPGmF7dgW0mTuS73oULi9cv1ELjbo6jnksA/DeAuUT0VwC3APh0pltVK4jwR3FBfX3Fwm+7j127iqeJi9q429PDz7Nv7+0Js8U1h0U9kvEDwJw57KTtOWClVos9SteOeOT3Xbu4R1JbW7DwL1rEccq+fc4yd9Qjr01L+Ht6HFH0q8JoX3iOPZbHQvhFPUGudWiI9/nQUHG/efeEJMaUJ/zujB8IF/45cxzhl4tcVOGfODG+43dfzEeN4u8typzF7iJ/UYgb9eSlVk/eoh5jzBMATgJwPIBPADjcGPN08KsahLSiHolahKiOf9s2R2QE98jUkSOjRz2A93y38jnl4ASKHb8t/BIZhTn+vXuB5593lvk5/iBnLZ8rLOMH+Dv6n//h35cs8X6uuNO2No4Rxo5N5vjt707KU9vbabtSGW8RV/iHh4svllOn8qNfz54g4Y8S9QDxHb/Ux7cRYY4SbYVFPV6EOf68RT3y/nmLeohoFIDLAFxsjFkNYDYRvSvzLasFxAmHCb8UYQuKeuyDNKrwy/vazxVHJgXKWlvjRT1eVTG9hN/P8csFRBoD3bf6gJPJStyzbx876LhRj9tJe2EL/y23cFfTBQu8nyvvf+ihfOfS2ZmsO2dc4U/i+GWayDhRT1sbXyBE+OXinpXjt2ffEuKUvAiLeryQc0xmOrPJY+OulG/JYdTzcwB7ARxX+LsbwLcy26JaYWDAOYGiVJAcHg6OeuI6/v37nYPFjoXc7QVRoh5b+A88kB9t4feKenbtck4sEX77zmH6dC4098EPlr7fIYfwa6Vnj7w+y6jnwQd5noPzzvN/rlykDj/c2R6vqCesgdIt/O47Ez/hjzMLmnz3EvV0dvJ+8BP+PXv4uJoyhY/Xffv4O+7qKnXGfsSpze9VpweIJ/xJoh5Zv4yJsfFy/BWKVgJpb89f1ANgrjHmOwD2AYAxph/BE7E3Bna3tjDHL0LgFn6/qCdKr54dOxyhsIVffpeTOSzqGRxkEZDnjxrFYiTCv3+/81ndjl8Kynk5foBnDLL/FkaM4O6S4vjdBdqEqMIfJeq5/np+3w98wP+5buEfN463bXiY3ytqd84kjt/+PFGQi7BcLInY9QdFPSL8AI9BidqVU+jq4u86ylSeaQh/kqhHzjGvi5m7cbelhe8Mqo0t/Dly/HuJqB1Or565AGIcoXWKxB8HHBBd+N2TP4dFPcPD/uu07zJsodm1i7NVOWHCoh73GAKgeKLzjRudhjh3xi/zs4ordgt/ENKzx5jSAm1CGo5/3DgWxe3bgbe/3b9Hj/3+tuPv7S1tCAxr3I0j/DKAC4gX97gdP8Axjl/dfLfwb97MUU8c4Z84kR+jtGml5fiTRj1ewt/ezhet4eF8zL4l2HciORL+rwP4HwAziegXAP4M4AuZblUtIMK/aFF41CPi6nb8+/axqHpFPUCwq7RPPrfjt+cgDYt63HcIQPF8txLzNDWVRj1ejj9qbLBoEbvOjRuTO/4ot+tNTY44BsU8AE868u1vA297G/8tUY+7LnzajbtJhN/t+AEuwx1V+Ddt4u84asMu4LQbRYl7yhV+OTfKiXrc2Ps5SftBVtifMQ9RDxE1ARgP4EzwZOu3AVhijLk/8y3LOyKIRx7JwhXUPc0v6gH4IPTq1WO/zgs/4XeLb1jU4x48BjjCb4xzgTv44NKoZ8oUvlX2i3qCkAbev/yltDKnkIbjBzju6ewE3hXSJ6GtjSdfkfVJ1OMW/vZ2p93GCz/hd3fnHDWKf9wzskXBy/HPmMGxnFdbgVv4n3mGxTeJ44/SwFuu8IdNwuJHmOMHKj7pSSj2Z8yD4zfGDAP4F2PMNmPMH4wxdxpjEs68XGesW8cn0dSpfKLZoujGS/jtg9DP8QcJv1/U476IRI163MI/MMBxgAj/woWOy9y/n99z7Finvg0QT/jf/GYu3/zpTzuNvEmjnrCT5QMfAL761fgnVWcnfzYRKdvxA/6uX76PtjZH+O08WV4v4wnScvwzZvB7e43eFeGXqEvKjOTV8SftdSPfUZjjz5Pw29uRB8df4B4i+jwRzSSiCfKT+ZblnXXruAeM1/R+brwyfnvmHa+M336dF3Ecf5SoxxZsu0vna6+x8501q7QXUUdHsfC7LzpBjBzplPmVKpZZOf5vfhO45JJo22XT2emUeZDtAcInY5Hv7cADHeF3n9xEpcIfp2xDby/HWPb+lujNq56OCH9HB7/fI4/w8npz/DKAMEj4KzyxeSjyGZubw4vlpUQU4f8ogE8BWA7g8cLPiiw3qiYQ4fea0NuNX8YPODGRu1cPkFz403D8AAu/fM7OTmd6OHuSDSllLPPtRhV+gOus/Nd/OZGJWyjSyPjLQS5E0lPGjnoAf/GS723WLEf47W0kcnpPAckd/7hxLP5CFOEn4jvVuH34geiOf/9+Pq6qIfwAH4Nex2Heo54Kbk/o5cUYM6cSG1JTSPb99rc7jj/IBQVl/DIZh1fUE1S2oafH6U/ujnpkrl0gXPj9Mn7AEf6ZM50subfXcfi243ePH4jKu97FDap33VXqdtJy/EmRGMUt/FGjnlmzeGIbL5EZPbo84bfLNQhRhB9g4X/1Vf7u3PFaEG1tvI4w4ZfjthpRD8CVZGXkuE1eHb9sR56En4jO9Fi8A8AqY8xmj//VP9u384Eb1fEHZfwi/HGjnm3buCtpd3ew448a9djv39nJblKinmOPdQSit7fY8Y8dy1GI1wUkKl/6Ev+4SSvjT4p8ZhkUFdXx79nDmf7UqXxcyKhZmx/+kBvMgeTC7xbtqVPZ0UcRfoAvTH7z2foRZfSuX4E2oDKO/73v9V7uzvi9tq8a2F2vK0SUQOkC8Kjd+wp/nwzgYQAHE9E3jDH/ldG25Rdp8Jw5M1rG7yWuURx/WNQzYQJfhNzdOe33SdKPH2BReO45dncS9QBOlUt5zdixPLF6UscfRFub0++6ySOVrFbUE8Xxjx7NIjk8zBcOt/Cfc47ze5JePXZlTqGlhc2A9DgT9u/nY8At/HEadoUo9XqChF8aud3Cv28fl/2WKSDLEX4/3FFP0JiOSlKFqCdKxj8M4DBjzFnGmLMALAAP4DoWwBez3LjcIsIvgkgUHvU0NRULVFrC7x7lGzfj37WLt999gs2aBTz8MP/ujnrcjr+cqCeIsBGt1Yp6ojTuivAD3Lc+6KROK+oBnC6d7u0BvB1/XMp1/NK+4Rb+228HTjkFeOUV/juLWjoa9bxBFOGfbYyxpuzBZgAHG2N6UCjj0HDYwi/FvMKiHmlYE4KiniiNu9u2OXVWxLXv388nVNxePaNHlzrqWbOKe6fYUY9Xxu9351AOYYJYKeGXqMc9ACuocdcWfq9J522S9OrxcvxAPOGvhuMHvGvyy8AzOR8q4fjzIvw5jXoeJKI7Afyq8PfZAJYT0WgAvVltWK5Zt46/JGmcc0/v58Zdix9Ip3F3woTi2v32ZBtClKjHK5e33eDMmcUzBckUjCL8/f1Ov/JKC39Li3cMlAYS9fhl/FEdf9j0kGk6/pkznbl+7e0BnONKIo5qOH7A2/HLeWBXQwXSFf68j9zNmeP/FLhC5yIARwG4GcCnjDF9xph/yHDb8sXPfw4cdxz3f37tNXZWIjgTJgSfDF7iGiT8I0dyDxc/xy+VOd1Rj5frjiL8XmItokDENXm8HL9EPYDj2Cop/O5ukmkzYgSLpfTjj9Odc/Rop+EfSFf4BwdZuPwc/86dxWWz3cK/eDF/v0uXRns/m4kT+RjYF3Czn0T45S7Cnv8AaKyoJ0+O3xhjiGgFgB3GmHsL9fnHANgV8tL6YetW4HOf4xGRxx/P4nb00c7/u7rCM36343dHPW7BDCrN3NvLXUol6pHZsrx61owc6cw369WDw6++jgj/1KnO3L2trfzeROy0W1urK/yDg9mfLJ2dzvcQp3F3/HjH8QPpCr/cXfpl/AB/H/LduIV/3rziiXbiYPdik8jITZjwjx4d7vgbMerJk+Mnoo8DuB3AtYVF0wH8LsNtyh9XXMHO+G9/A84/ny8A8+Y5/08j6nGLb5Dwy3u5ox6vBlYRRr+c3y/qkfzXHhMwfrzTq0dOarfwJ+nO6UcU4c/6ZBFXLRNmANEd/6hR0RruRo7ki2lU4ZcL/bRppf/z6ssvx5F7NqwkRBm9K8LvdyxEcfxZNu7u2VOZYycqeezHD456jgHwCAAYY14iopz0g6oAzz8P/OQnwIUXctRz3HHARRcVC2JY1NPXV+qE5UveupVPfHHVQlTh94p63N05AX93vHu3d7c2mRJQJmYBnDLFra3O55EcfP16FsY0HXgeHL98Ptt5Ru3OScRCGda4SxRvFi65yEpZbJsg4XebjyREGb27cSObBL9a96NGOfGZIAbInvgGyCbjl/fIi/DntHF30BizlwoxARGNQKE2f0Nw6aV8oF5+ubNs8eLi53R1OaUXvA723budOVEFOeiGh73jEbu3jhsRfnevHi/HLxeUIMd/0EGly4mAr3ylONISxz96tLfjt8tBp0G1M37Acfy2AElf9DDhBxzhD9tOuyZ7GEHCL3cBWQl/FMf/xBNcdtuPqBl/2hOlNDXx+VDhSU9CyWPUA+ABIvoygHYiOhXcu+eObDcrJ7z8MnDnncAXvhA82EMGcflV6PSKeuQgBLxviYMcv5x0dtQj8/q612c7fi+Cauh/5StclkIQx79zp3NxcQt/muTB8Yvwu0/KoMlY3MLv9Xo3fo5/eBh49tniZevXsyB6ZeytrXysVsvx793LE+z4TWgP8PFmNz7v2+ecO15zHKdJW1vFJzYPJaf9+C8DsAXAKgCfAHAXgH8LexER3UhEm4lotbXsu0T0PBE9TUS/JaLOhNtdGcRZS+14P+Rk8HNBXsIPOF+0l2C6hf/SS4Ff/rJ4uyTq2b+fT7igjN9P+P0yfi/Gj3cGcLkdf9wCbVHIQ8bvFfUAzmQs+/bxBfLll3n58LAzvy1QvvD/5jc8TaUMbAJY+A84wN8Nu/vyywUqTeH3O9afeYaPBfddsc3UqdxFVorz2euyHX+aMY9QhWkOQ6lC1BMq/IWa/L8DcJEx5mxjzPXGRJoV+iYAp7uW3QPgCGPMmwC8CMCjQEuOkBMmrFEsrGyDDJJyE1X4jQF+/GOu72K/z/jxjmjv3u3t+IOiHrlLiCrYnZ3slrwcv9/nKIe8Rj3y9549wI9+xEXmfvMbXi5xTVrC/9JL/D09/7yzbP1675hHmDGjuGxDmo5fJo/xc/wrCoV7gxz/9Ok8FkRyfsn3geyFv63NeY8s1p+EPEU9xFxORFsBPA/gBSLaQkRfi7JiY8xyAD2uZX8yxhRG/+BhADMSbndliNrAFCT84gCDZgQKi3p27eJ1PPYYu9xt21iQmpuLB3t5decMcvwDA3y3ENXx21GPCL49IjnNHj1AvqIeL8f/6qvA17/Of4trdYtsVOG352C2kTxf7ihkWZjwu6OeESNKOxAkJWj07ooVvM+82o0E2Xb3aN0JE7KPetrbNepBsOO/GMBbACw1xnQZYyaA6/O8hYg+l8J7fxTA3X7/JKILiWgFEa3YYjuCShLV8Qfd/rodoE2Q47d760j3vYEBbjiTUbvyPICfu3s3b6sdAQR15/S6Qwhi/Hi+UGze7GwzkXMRqMeM3y/qaW8Hli93Yh05Rt1dJ8t1/CKO7qgnSPhnzmRxk23xixqTEjR69/HHOeYJauR39zySi8j8+cW9erJ2/HkR/pxFPecBONcY88YRZ4x5GcCHC/9LDBF9BcAQgF/4PccYc50xZokxZskkKY1QaUS0y4l6gm6zw6IeEWYpEgZwBcOeHudiY0c9Xg214vK8HH9c4Rf3OzxcHPFUU/gr1Y/fS/gB4LOfBebOdcQrqeP369Uj4iiOX6ZWDHP8gHPRSFv4/Rz/4CDw9NPBMQ/g7/jnzatM1JM3x5+nqAdAi9f8usaYLQBakr4hEX0EwLsAfChiW0H1cM+16se4cf4VOr1m3xJkvX7C39/PImtXhxThdzt+EX73uoKinriF1eyRovZrqiX81cz4J0zg3jNf/SqLexrCHyXqkb+9Bm8JbkddKce/ahU3dgc17AK835qbnc8i+27uXN7WffuyjXqk3ERehF/Om7TPnwCC+vEHlHQM/J8vRHQ6uJTzScaYkJkYMmTTJj7AgnJIILrjb25mURTHv2uXU/HSa75dQQ48v4wf4IuPRD3veAdw//18oZk7t/h5EvW41xUU9cSdPMWuDWMfpBKHpH3gyrbnIepxi8QPfsACMm4cD3Z78kle7hZ+GegnF2o/vIR/aIiPVYCjHmOcYyEs6gG4DUK2qRKO//HH+THM8Tc3c88e2/GPH+/cxe7YwedeFhOl2N9jXoR/0iSegW7Zsoq9ZZDjP5KIdnr87AKwMGzFRHQbgIcAHEJE3UR0AYAfAugAT+C+koh+msqniMsHP8izH116aXDp46iOH3DKNjz1FJ+UV1/Ny8uJeuT1Gzfyxecd72Cn9fLL0R1/mlGP7fgrEfUQsbDnsXF31iynbEeQ4z/4YG7wPN3dwc2Fl/BLl8dDDuEG9Z6e4MFbwoEHsumQdoEsHH9vr1OlVVixgo/LKOWep08vdvyTJhUXAhwYyC7q8fq92rz97el3jgjA1/EbY8oaMmeMOddj8Q3lrDMVdu8GHnyQT9zvfQ/41a+Ae+7xnqMzToXACRO4y927380C/NxzvDxI+MOiHnn9xo3skGxH4JXx795dXBESiBb1lOv4RfizOHCDShlUM+O3mTjRKVft9X2HRR+yfvfnFGFctgx44QW+4EcR/pYWdv228MudSxpMmcJ3H5s2FW/HihXhDbvC9OnOwLQtW3gfyr4Wx59V1CPkSfgrTEaFzHPM8uV8i37ttfz7unXAbbd5P3fPHj44otR77+ri2/1t24pvY4My/qCoxxb0DRt4nfPmOSOIxfHbF4ggxx/UqydOP36hEo4fCBb+SmT8fr16bCZOZCGU+XWB+A7b63PKMXTCCfz4yivOCOmwfT1nTnGDcJqOXyq3SpQE8LavXh0e8whhjj/Lxl2v3xuMxhP+e+/lL/yEE9hJzZzJg2S86O+PXtFw0iR2Orfeyge/ZLFRMv4ojn/aNF6/iIBb+MXx+2X8Xo4/bsZvu8ZKNO4CwaUMhoayF/62NuCMM4LzV+l1tnVrecLv7tVjO37AcfxBbl846KDsoh4v4X/2Wf4+jjoq2jpmzODjb9euUsefZdSjjh9AIwr/PfeweMqXPm8e8Pe/ez93z57oB99XvgLcfTeLxPTppcJfTtQjjh9whF8iHamLX06vnqjC39xcWqrB/r2Swi+fJ+uTlwj43e+Ad77T/znSc6dc4R8Y4DsHYf16/n5nz+aLSxzhnzOHDUN/f2WEf80afjz44GjrsLt0iuMXYyGOP6taPV6/NxiNJfyvv863o299q7Ns/vx0HP/BBwNvexv/Pm0aH8yDg9Ead4N69WzaVFzd8z3v4dotb3pT8XPD+vH7RT1NTfGclTTwVtvxy7IKDnrxRYR/yxb+vqXMchza2py7GGH9ej6WmprYwcd1/ABPuJK28I8Zw3ecXsIf1lNOkM/w7LMcvdqOf/v2yjTu5uHYqRK1L/yXXAL8+tfRnitzkZ56qrNs3jzO5WVQh00cx28jfaw3boyW8Qc5frkbkXXOncv9paWvNsAn4tat3iWew6KeMWPilVKWkzMvjj8PJ68d9chI3rjlqb3GLNgiP2cOi+vGjdEdP8CvsYvGpcWsWaXCP3ly9GNAPsPKlfw4aZJT1lu6sGbZuNvWlm4J8RqjtoX/0Ue52+Ttt0d7/j33sFOxa4VLbx6vuGfPnmSzFslBvWEDu62WFu86KVGEX+5G3PX8bcaMcQZ5xe3OGbcnTmenEy8Jp5wCnHuud8+ocqkF4ZfYTRx/EpG1pwUUbOE/6CB270ND8Rz/s89yfJTG7Fs2XsIf1e0Dzmd46il+nDiR72zGjSud3D5NqlAXJ4/UtvBfdRU/RqnlYww37P7jPxb30hGx8op7kvYsEHcuwu8nBEEZvwiyXJCChH/0aOdkcQt5czP/+EU9cYV//PjSgTVz5nCjdlaZbDUz/ijIjGSS8ScRfnsicICPV7fwC1GEf/JkFvvVharoWTl+aZN4+WVnUGEURo1iE2E7foCX2SPV06YK5RHySJQZuPLJ2rWO0w+aBk544QU+keyYB+ATisjf8QcJrh8i/OvX+5dkBoCzzmJBdve9B0odf9AQ/TFj+PMB3heR1lZ/xx83njnkkOLaQVlTCxk/wMKVhvDL59q5k9dlRz1CFOEn4tesWsV/ZyH8fX3chbWjg8tAxxF+gONKuTBJO0lnp2NismzcbXDhr13Hf8017NxPOy2a43/wQX48+eTi5W1t/l06kzr+ri6OWMIc/5w5wJe/7J01trTwT08Pi5vdh96NPaORl4P3E/6g2bf8+Na3gAceiPeacqiFqAdg4Son6nELv3ugVlzHD/DxJQMJsxB+gF3/2rXcvhRX+O3PIY5/3LhsHb8KP4BaFf7eXuBnPwPe/37O67duLe4G58Xq1Xzwex2cfl06k2b8ROzQRfiTjmiVk3Xq1OCGKPuk9nLwI0d6Rz29vfEdf9oTqodRS8KfpuN3C/+MGU5s5zXlohcHHeTspyyFXwaKxcn4AeeztbY629fZ6dS80qgnM2pT+G+5hWOKf/1Xdgr2tIN+rFoFHH649yhcvy6d5YwenDaNT95yutLJ64JiHqD4whLV8Q8OshtcsCDZtlWKWsj4gfSjHrfwjxjBYhs05aIbOx7KSvjXrnW6ciZ1/DL4ESi+s9WoJzNqU/iff5575xx1VPHgmSBWrwYW+tSWmz/fu0tnUscPOIO4gjL+METEw9oZbLGPmvE//TT3n166NNm2VYpayfjLjXrcvXqkpLIdhyxcGH2AFFDswNMW/q4uZxayNWv49wMOiLcO+WxyDgPFwq+OPzNqU/jtiUgkGwzK+Tdt4v8fcYT3/6XKoh33GFO+4w/L+MOwo54ozwO8Hb9X1PPYY/xYq8Kfx6inv5+PszR69axfz8e4LVA33gj8v/8XfZ1ZOn4iHlEswi+dJOJgO34ha+FXxw+gVoV/2zZH+KM4fuk5EOT4geK4Z98+nmYwqeOfNo3jp02bys/440Q9UR3/o49yl78DD0y2bZWirY0vWsPDxctffJEf0+6fnhQRrx070ot63I24EyYUi2QYWQo/4HTpjNuVU5BBiLbjt+tBadSTGbUv/FEcv3Rp83P8Xl06o07C4oectD092Tt+EX73wCrBS/gfe4zdft5HL8oJam//Qw8BX/sa17iXu7VqY4tXVsIflzFjnPMjK+Ffu5aFP27DLlAdx69RD4BaFn6pThnV8U+c6JQ0duPVpTPOJCxe2C69UsLvd2fhjnpkvoC8xzxAqSBu3MjjH2bO5EFjeblw5VH4Acf1ZyX827fzuZLE8U+cCBxzDHDccc4yjXoqQm0O4LIz/o4OFrYgxy8Nu0EiMW9esfCX6/ht4c866pHn+XXNbG3lCEJ44gluw6hF4T/3XP4sf/xj8Yxg1cZ2reU27vb1AZs3pxPDHXQQ391lIaLSswdIJvxEwCOPFC/LuleP7Ics9kcNUXuOf+/e4pmmiIqnvnMzPMzC7xfzCPPnF0c9eXD8cXv1+F1g3FFPrTTsAsXCv2sXDx77whf822uqRZqO/+mn+cJsV2BNyrJl4aYnKeUKvxfq+CtC7Qn/tm38aJc5kD7UXrz6KjuoMKGYO5fXLSNgy3X89ixJSYV/8mSui+NV0sFGBN/P8bujnscec2q85x1bEKVve15yfZvOTqd/fZLvW9pm+vudidujTmoSxEUXOYXQ0kaEv6mp+CJQDrbwZ9FjS4UfQC0Lv2T8gNOH2ouwhl17HYDTl79cxw84rj+p8H/uc8DDD4e7NVl/HMdfC24fKBZ+6dtul6TOC01NzgU6yffd1MQX6IEBFv4JE7gdI89MnerM7+tVfTYJ0qsnq7LJLS1c5uXYY9Nfdw1Re8Ivw7mjOn7pynn44cHrtWf/Acp3/IAj/Ekz/rFjgcMOC39emOO3hX/LFp6SrxaFXxx/HoUfcMxD0gu9jFlYuZLdfl4arv1oauJ2iLRiHsCp/JpVBk/E7UNnnJHN+muE2mvc9Yp6whz/rFmlpYTd2PN9Auk4fumVkUWPCps4vXpWrODHWhR+cfxhjd3VIg3h37WLj9lPfzq97cqSa64JLiAYlxEj2MA0eONr1tSP8Pf28qCrlpbi50dp2AVKhT9Nx5+18Efp1SOOX+qfH310ttuUFm7h7+rKryiU22e+vZ2/n8HBdPL9ShA0F3FSZLIfJTNqL+rxyvjlhJP/Cbt2cV2fI48MX68Iv3R7zEPGH5WwjH/kSEf4//53rqkSdgeUF9zCn9eYB0jH8UtDbK0IfxZ0dub34l4n1J7w9/SwkNknl98grnvv5anqTjstfL1ZZPzHH8+9Z9Lq8eDHiBHAF78IvO993v9vbXWinjVr8tkrxg93xl/vwr9/P4tenGJs9UZnZ8P3usma2ox6urqKG778yjbceScfRMcfH75et/Cn4fiXLuWG1Epw5ZX+/2tt5fEMQ0Ps+N2zkOUZt+PPc9vEkUdyrfykmbd81je9KXrp5XrkYx9zzj8lE2pX+G28HP/wMPCHP3A9lyh54YgRHJW4hb8enId0tevtZddci45/xw6+sOfZ8Z9zDv8kRT5rI8c8AHDeedXegrqn9qIeL+H3cvyPP86VMd/1rujr7ux0Mn4pyZz3LnVRkIEwMg1fLQq/zPKUZ+EvFxV+pULUnvD39BQ37ALOhcB2/Hfeyf2MTz89+rrHjSt2/Hkp+VsuIvzPPMOPafa7zhoRQymnkUbhsrwisaIKv5IxtSf8Xo6/pYXduu3477iDs/2wcgc2nZ3Fjbv10rNAop5nn+XHWhJ+uWiJ8Ne7429uzl8dIqXuqK2M3xhv4QeKC7WtX8/D3oMaPL3o7OSyv0B9Ov5nn+V9l6eqlmEQ8favXct/17Pwn3QSf9Z6aFdSck1tOf7du3mQlpfwT5rkOP4//IEf4+T7AEc97oy/HrCjnlrK94W2Nv7eOzpqZ/xBEj75SeCmm6q9FUoDUFvCL3V63Bk/4Dh+Y4Cbb+Y4Y8GCeOu3o556cvwS9bz+eu0KP1Df+b6iVJDaEn6vcg2COP4HHwT+9jeubBm3R44Iv0y0Xi/Cb5e3rWXhr+eYR1EqSGbCT0Q3EtFmIlptLTuHiJ4homEiWhJ7pUHCL47/29/mOvYf/Wj8je7s5JGTfX3s+Ost6gFqq2FXUOFXlFTJ0vHfBMDdl3I1gDMBLE+0xjDHv3cvl1y9+OJkom2P3q0nx2/XSq9lx69Rj6KkQma9eowxy4lotmvZcwBAcSOY11/nx7CMH+DGv4suird+wS7UVq+Ov5aFXx2/oqRCbrtzEtGFAC4EgEUjCpvpVZlTEOH/1Kcc5x4XuzRzPTl+Ef6xY4vnhq0VVPgVJVVy27hrjLnOGLPEGLOkeWgIWLeOhb+jw3uatxNP5Ijn859P/qa28NeT45f9NW9ebZag0KhHUVIlt46/hHvu8R+8BbCbvfrq8t6jXjN+cfy1GPMA6vgVJWVqQ/hbWoA//YkHcHnFPGkhjn/zZq7uWS+Ovx6Ef+TI2oypFCWHZNmd8zYADwE4hIi6iegCInofEXUDOA7AH4joj5FWNnYsT6qyZUu82jtxEccvZRvqxfF3dQEnnBCvYF2eOOig2ph8XFFqhCx79Zzr86/fxl7Z2LE8ocnOncDZZ5e3YUG0tfHPhg38d704/pEjeWBbrfKtb/H4CkVRUqE2oh6ZRHxoKFvHD7DrrzfHX+s0NfGPoiipUBtnU0uLM2F61sJvV+isF8evKIpiURvCDzgTpmfZuAsUC786fkVR6pDaE36ZZjEr7Aqd6vgVRalDakf4TzkFuO464Iwzsn0fe9SvOn5FUeqQ2mjcBbhx7+Mfz/59pC8/oI5fUZS6pHYcf6WwhV8dv6IodYgKvxt1/Iqi1Dkq/G4041cUpc5R4XejUY+iKHWOCr8bEX4i7/LPiqIoNY4KvxsR/lGjtCiYoih1iQq/G8n4tWFXUZQ6RYXfje34FUVR6hAVfjci/Or4FUWpU1T43YwaBTQ3q+NXFKVuUeF3Q8SuXx2/oih1igq/F52d6vgVRalbaqdIWyU57DBg+vRqb4WiKEomqPB78bvfaR9+RVHqFhV+L5qbq70FiqIomaEZv6IoSoOhwq8oitJgqPAriqI0GCr8iqIoDYYKv6IoSoOhwq8oitJgqPAriqI0GCr8iqIoDYYKv6IoSoOhwq8oitJgqPAriqI0GCr8iqIoDYYKv6IoSoOhwq8oitJgZCb8RHQjEW0motXWsglEdA8RvVR4HJ/V+yuKoijeZOn4bwJwumvZZQD+bIyZD+DPhb8VRVGUCpKZ8BtjlgPocS0+A8DNhd9vBvDerN5fURRF8abSM3BNMcZsBABjzEYimuz3RCK6EMCFhT93E9ELEdY/EcDWBNs1DsCOOnqN7gcmyX5I8j5JX6f7obKvqdR+yNO+m+W51BiT2Q+A2QBWW3/3uv6/PeX3W5HwddfV2Wt0PyTcD0neR/eD7oda2nfGmIr36tlERFMBoPC4ucLv78cddfaapOT5M1VqPyR9H90PyV9Xb/sh7/sOVLhqZAIRzQZwpzHmiMLf3wWwzRhzJRFdBmCCMeYLKb7fCmPMkrTWV6vofmB0PzC6HxjdDw5Zdue8DcBDAA4hom4iugDAlQBOJaKXAJxa+DtNrkt5fbWK7gdG9wOj+4HR/VAgU8evKIqi5A8duasoitJgqPAriqI0GLkXfp/SD4uI6GEiWklEK4jomMLykUT0cyJaRURPEdHJ1mvuJ6IXCq9ZGTSGII/47Icjieihwue9g4jGFpafSkSPF5Y/TkSnWK9ZXFj+dyL6TyKianyepMTZD9b/DySi3UT0eWtZIx0PH7I+50oiGiaiRYX/1ezxEHMf1K02JCJJH9BK/gA4EcDRKB4P8CcAby/8/g4A9xd+/xSAnxd+nwzgcQBNhb/vB7Ck2p8n5f3wGICTCr9/FMA3C78fBWBa4fcjAKy3XvMogOMAEIC7ZT/Wyk+c/WD9/9cAfgXg89ayhjkeXK9bCODlejgeYp4TdasNSX5y7/iNd+kHA0Bc3TgAGwq/LwDXAIIxZjOAXgB10X3LZz8cAmB54fd7AJxVeO6TxhjZJ88AaCOi1sLYibHGmIcMH/G3oMbKZsTZDwBARO8F8DJ4P9QNcfeDxbkAbgPeGEtTs8dDzH1Qt9qQhNwLvw8XA/guEb0G4HsAvlRY/hSAM4hoBBHNAbAYwEzrdT8v3Mp9tZZuaQNYDeA9hd/PQfFnFc4C8KQxZhDAdADd1v+6C8tqHc/9QESjAXwRwBU+r2vE4+H9KAg/6vN48NsHjaYNgdSq8P8fAJ8zxswE8DkANxSW3wg+eFcA+D6AvwEYKvzvQ8aYhQCWFX7+uZIbnBEfBfApInocQAeAvfY/iehwAP8B4BOyyGMd9dCf128/XAHgamPMbo/XNOLxcCyAPcYYycTr8Xjw2weNpg2BVLpIW1p8BMBnC7//CsDPAMAYMwS+EAAAiOhvAF4q/G994XEXEd0K4BjwrW3NYox5HsBpAEBEBwN4p/yPiGYA+C2A84wxawqLuwHMsFYxA05MVrME7IdjAZxNRN8B0AlgmIgGjDE/bLTjocAH4Lh9oA6PB7990GjaEEatOv4NAE4q/H4KCl8gEY0q3N6DiE4FMGSMebZwezexsLwFwLvAt4Q1jfQ+IKImAP8G4KeFvzsB/AHAl4wxf5XnG66MuouI3ly4nT0PwO8rvd1p47cfjDHLjDGzjTGzwS7v28aYHzba8WAtOwfA/yfL6vF4CDgnGkobwsi94ycu/XAygIlE1A3g6wA+DuAaIhoBYABO+ebJAP5IRMMA1sO5ZWstLG8B0AzgXgDXV+xDpIDPfhhDRJ8qPOU3AH5e+P1fAMwD8FUi+mph2WmFRq3/A54kpx3ci+PuinyAlIi5H/xotOMB4B4w3caYl12rqtnjIeY+qFttSIKWbFAURWkwajXqURRFURKiwq8oitJgqPAriqI0GCr8iqIoDYYKv6IoSoOhwq8oLoiok4guKvw+jYhur/Y2KUqaaHdORXFBrrmiFaXeyP0ALkWpAlcCmEtEK8Gjwg8zxhxBROeDq1c2g8tdXwVgJHgw0CCAdxhjeohoLoAfAZgEYA+AjxdKCShKLtCoR1FKuQzAGmPMIgCXuv53BIAPguu5/Du46NlRAB4ClzwAeFLvTxtjFgP4PIAfV2KjFSUq6vgVJR73GWN2gWvc7ABwR2H5KgBvIqIxAI4H8Curum9r5TdTUfxR4VeUeAxavw9bfw+Dz6cmAL2FuwVFySUa9ShKKbvAtdxjY4zZCeAVIjoHAIg5Ms2NU5RyUeFXFBfGmG0A/lqYxPu7CVbxIQAXENFT4Ckfz0hz+xSlXLQ7p6IoSoOhjl9RFKXBUOFXFEVpMFT4FUVRGgwVfkVRlAZDhV9RFKXBUOFXFEVpMFT4FUVRGoz/HyC6S/RFiieDAAAAAElFTkSuQmCC\n", + "text/plain": [ + "<Figure size 432x288 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX4AAAEHCAYAAACp9y31AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABL70lEQVR4nO29eZQjd3nv/X1KpSptvXfP0rPbxOt4ZWyzr3HeG+MAAXOvCQnkwglJ4OZeuDcL4Q1Ltnu4CSRvOLk5xCFASIAcMA5mCUnAwZjFxhmv4xUzHs9MT89M791aWiqV9Hv/qPpVlaQqqaSWWtWt53POnOmulkqlavW3nvo+z+95SAgBhmEYZnBQ+n0ADMMwzObCws8wDDNgsPAzDMMMGCz8DMMwAwYLP8MwzIDBws8wDDNgqL3aMRF9CsDNAOaEEIftbW8C8GEAlwK4XghxNMy+JicnxcGDB3t0pAzDMNuTBx54YEEIMVW/vWfCD+AzAP4SwGc92x4D8AYAf93Ojg4ePIijR0NdIxiGYRgbIjrpt71nwi+EuIeIDtZte9I+mF69LMMwDNMC9vgZhmEGjMgKPxG9k4iOEtHR+fn5fh8OwzDMtqGXHv+GEELcBuA2ADhy5EhDQ6FyuYyZmRkUi8VNP7bNIpFIYO/evYjH4/0+FIZhthGRFf5WzMzMYGhoCAcPHtyWOQMhBBYXFzEzM4NDhw71+3AYhtlG9MzqIaIvALgXwMVENENE7yCinyeiGQAvBPANIvrXTvdfLBYxMTGxLUUfsBLgExMT2/qOhmGY/tDLqp43B/zon7r1GttV9CXb/f0xDNMfIpvcZRiGYWqZWS7gO0/NbXg/LPwdsrKygr/6q7/q92EwDDNAfPbek/jVf3gAGx2gxcLfISz8DMNsNvmSCcOsIlsyN7QfFv4Oed/73ofjx4/j6quvxnvf+168+tWvxrXXXosrrrgCd955JwDgd37nd2ouDh/+8IfxsY99DNVqFe9617tw+eWX4+abb8ZNN92E22+/vV9vhWGYLUKxXAUALOWMDe1ny5Zzevn9rz2OJ2bXurrPy6aH8aGfuzzw5x/5yEfw2GOP4eGHH4ZpmigUChgeHsbCwgJe8IIX4LWvfS1uvfVWvOc978G73vUuAMAXv/hF/Mu//AvuuOMOPPfcczh27Bjm5uZw6aWX4u1vf3tXj59hmO1HyawAABbzBg5Opjvez7YQ/n4jhMD73/9+3HPPPVAUBWfOnMH58+dxzTXXYG5uDrOzs5ifn8fY2Bj279+PP/uzP8Ob3vQmKIqCXbt24ZWvfGW/3wLDMFsAGfEv55tH/EII3PTx7wf+fFsIf7PIfDP43Oc+h/n5eTzwwAOIx+M4ePCgU39/yy234Pbbb8e5c+dw6623AsCGEzMMwwwmMuJfaiH8xXIVT54NdkHY4++QoaEhZLNZAMDq6ip27NiBeDyO73znOzh50u2Eeuutt+If//Efcfvtt+OWW24BALzkJS/Bl7/8ZVSrVZw/fx533313P94CwzBbjGLZtXqasbpebvrzbRHx94OJiQm8+MUvxuHDh3HdddfhqaeewpEjR3D11VfjkksucR53+eWXI5vNYs+ePdi9ezcA4I1vfCPuuusuHD58GBdddBFuuOEGjIyM9OutMAyzRSiZdnI3X2r6uLUiC3/P+PznPx/qcceOHav5XlEUfPSjH0Umk8Hi4iKuv/56XHHFFb04RIZhthEy4l/KNxd2jvgjys0334yVlRUYhoEPfOAD2LVrV78PiWGYiOOUc7aK+Fn4own7+gzDtEvY5G6riH9LJ3e3e3XMdn9/DMO0hxPxFwZU+BOJBBYXF7etOMp+/IlEot+HwjBMRHA8/hYrd9fWm7d02LJWz969ezEzM4PtPJZRTuBiGIYRQqBkVhGPEfJGBcVyBYl4zPexq+tlZPRged+ywh+Px3kyFcMwA4Ms5dw1ksDppXUs5Q1MjyZ9H7tWLGM4ESzvW9bqYRiGGSRKtr8/PWKJfbME7+p6GcPJ4FndLPwMwzBbgKJd0SOj/GbCv8bCzzAMs/WREf/uEavgo1XEP8LCzzAMs7VpJ+LPFk0MJ1j4GYZh+oJZqeLNt92Hbx47u6H9yFLOHUM6YgpxxM8wDBNVvveTBdz77CLueOjMhvYjq3pSmoqxVDywQ6dZqSJXMjGc5KoehmGYvnCnLfj/8dwSqtXOF5zKiF+PKxhPa4H9erJFa/EWR/wMwzB9IF8y8a+Pn8dkRsdKoYynz2c73pds15BQYxhLaVgO6NAp2zWwx88wDNMHvvXEeayXK/h/X2PN6PjRs4sd70tG/Im4gomMhsWAiF/24ueIn2EYpg985eEzmB5J4HVX7cGe0SR+dGKp431Jj19XY7bV4+/xOxE/Cz/DMExrDLOKX/ib+/B7XzmGEwv5De1rIVfC955ZwOuu2QNFIdxwwTjuP7HUcWNJb8Q/ntKwsl5GxSdnIBu0ccTPMAwTgvNrRfzw+CL+4b5TeNXH7sZ//8JDHQv1Nx49i0pV4PVX7wEA3HBoHIt5Az+Zy3W0Pze5a0X8QgArPu2ZZcTPws8wDBMCKZr/++evwE1X7MZXH5nFSqF5b/sg7j+xhH3jSVy8awgAcMOhCQDAfR3aPa7Vo2A8owPwX8QlPX4u52QYhgmBFM1Dk2m8/KIpAECu1Ly3fRALuRJ2D7vdMw9MpLBzWO84wVsqV0BkCf9EWgMA31r+1fUy4jFCMqBlM8DCzzAM4yD98eGkiiG7n33e6Ez4F/MGJjKa8z0R4YZDE/hRhz5/0axCVxUQEcZS1n6X/SL+9TKGE3EQUeC+WPgZhmFsHJskEUfaFv5csTPhX6oTfgC4/tA45rMlzCyvt70/7+AVud+giL+Zvw+w8DMMwziseUohHeHvwOoxK1UsFwyMp/Wa7c/bkQGAjiqGSmUr4gfgRPx+Hv/qehlDLPwMw2wnVgtl5Dv03VuRLZogAoZ0FUOJzoV/uVCGEMBkXcR/YCIFADi5VGh7n0XTjfg1VcGQrmIh17iIa61ocsTPMMz2waxU8fN/9QN88M7He7L/taI1q1ZRyIn4O7nIyFW1E3UR/86hBDRVwanF9iP+YrmChOombJ+3M4Mnz641PM7y+JtP1WXhZxhmy/DNx87h2YU85rLFnux/bd3tYy+HlWc78PgXc5YFU+/xKwph/3gKpzqI+EtmFXrclexr94/hkZlVGHaZp2SNPX6GYbYLQgh84rvHAQAFo9KT11grlh2LJ61Z0XW+1P5rSQum3uoBgAPjKZxc7MDqqYv4r90/BsOs4glP1C+EaDlvF2DhZxhmi/CDnyzi8dk1aDGlZx6/d1atGlOQjMeQK7W/gMuJ+OusHgDYZ0f87ZZ0Fst1Ef+BUQDAgyeXnW3r5QrMquCIn2GY7cEnvnscU0M6Xn3pDqyXW0fhz87n8Pw//BZOt2GrrNWNLEzrKnIdRPyL+RJiCvkK8IGJFApGBQu5xoqc00sFfPJ7z/rus1iuQPdE/LtHkpgeSeDBU67wh2nXAPRQ+InoU0Q0R0SPebaNE9G3iOgZ+/+xXr0+wzDbh8fOrOL7P1nA2198CKOpeCir5ydzOSzmDTx9LnwPfCvidxOjQwm1s+RuzsB4WoOiNC6ikpU9p5YaE7xff/Qs/ugbT/qWaRpmFYl4rWRfc2AMD51a8Ry/vQCtSS9+oLcR/2cA/Ke6be8DcJcQ4qcA3GV/zzAM05S7n54DAPzC9fuRjKtYDyH88uIQ1Lfej2yxXBfxxzoq51zMG05bhXr2j6cBwDfBW7IHqvtN1/Iu4JJcu38MZ1bWcX7NSnb3PeIXQtwDoL4b0esA/J399d8BeH2vXp9hmO3DXLaEkWQcI6k4UloMecNs6ZFL4fezVPyoVgWyJbMmMZrW1M6EP1dqqOiR7BtPggi+CV7ZiG3R55iLPhH/tftHAbg+v9uLP1rlnDuFEGcBwP5/xya/PsMwW5C5tRKmhqxEaUqPQQhXJIMo2D125rPhIv6cYUII1NTADyXUjlo2WBF/Y2IXsAap7B5O4JSf8NvjFZd92i2X6jx+ALh8egSaqjg+/1q/I/6NQkTvJKKjRHR0fn6+34fDMEwfmcsWsUMKv213eH3+1fUyflw3z9a1esJF/Gs+s2rTutpRk7bFXGOfHi/7J1K+q3eNSvAx+0X8mqrgij0jeKA+4u+jx+/HeSLaDQD2/3NBDxRC3CaEOCKEODI1NbVpB8gwTPSYz5Vc4desiLzgEeRPfu9ZvOkT99Y8x7F6Qkb83s6ckozefsRfLFeQK5mYzPhH/AACF3HJiH+pzuopV6qoVEVNHb/k2v2jeOzMGkpmxWkyNxSxlbtfBfA2++u3Abhzk1+fYZgthhACc2sl7BhOAACSWmPEv5ArYXW9DLPi2j/ywhA2uevtzCnJ6O17/DJaD0ruAsCBiTTms6Waixfg8fjrIn53+lajZF93cBxGpYo333YffviTRWR0FWqsubT3spzzCwDuBXAxEc0Q0TsAfATAjUT0DIAb7e8ZhmECWSuaKJlVTNkRdFpvFH7ZViHv2dZuclfuw5vczegqSmYV5UrzfIKXRXvV7kSLiB9orOyR7RfqyznlBaG+qgcAbrxsJz5482VYyhu4/7kljKWb2zwA0Px+YAMIId4c8KNX9+o1GYbZfszbfXl2DFtCmow3Wj2y1r5guJ0p5c+XCwbMSrVlFBzk8cv9j6aCI3gvQX16vDhdOhcLuGTXsLNdlnPWJ3edQes+Vg8R4e0vOYRfftFB/OD4QtPJW5KeCT/DMEw3mFuzIminqse2ery1/LKfjnexlYz4hQCWCgZ2DCWavo6fP55JuI3awgq/06cnoKoHAA7IWv66yp6gcs6i7f37WT0SRSG89KfC5UMjW9XDMAwDWIldAI5wS+H32jpZW/C97RUKnq/96uLrkcndGuGvG784ly3iuj/+Nh45vRK4H2nTNIv4R1JxjCTjbVg9tsfvE/F3Ags/wzCRpiHit8V43c/q8Ub8ZdOZm+s3sKSetWIZaS1WYwnVj1/8yfkc5rMlfOfpwIJELOYN6KriXKCCOOBT0lnyCL93gZqM+OvLOTuFhZ9hmEgzly1CVxVnYZVfHX/OifhrrZ59dhI1lPD7tDPO1I1flHcfzSL+hVwJkxm96bBzABhPa1gp+Ef2RqVac0dTKnPEzzDMADGXLWHHsCukfuWcUpi9i60KpYpTPRPG6snWdeYEPFaPbRvJVcCPzqwGtoxotXhLklBjTtJW4l2N7K3ld6t6OOJnGGYAmM+WahKzuqpAIbdqxzCrjjfuHZpSMEzsHNahxRQnUm/GWrHc0OMm48zdtRK/cj+LeQNnVtZ997OYLzWt4Zck4opj4UgMs+os/PKuP3CqekJU7ISBhZ9hmEgzly05NfyAVb6Y1lQn4vdW8tRX9aR0FRMZLVxyt64zJwBkNCn8bsQvHZxHTq/67seK+IMreiSJuH/EPz1qXeS8Cd6iycLPMMwAMbdWdGr4JUkt5pRz5nyE3zCrMKsCaS2GyYwe0uM3G1odyMViMrm7kDNw6a5haDEFj86sNOxDCBHe6vET/nIFu+wVyt7Vu045p9odyeY6foZhIkuxXMFa0XT69EhSWsyJ+GuE394mbaCk1mbEX5fcVWMKEnHFyR3MZ0vYM5pAPEZ4xEf4syUTRqXatIZfoseVhg6jRqWK3SOW8C97hL/EVg/DMIOCTKbWL75Kaqoj7n5Wj7wopEJG/EII3+QuAGT0uNPOYT5rtYe+at8ojs2solKtTfCGWbUrSagxlMyqkySuVAXKFYGxtAZNVeqsHk7uMgwzIMxla2v4JWlPxJ/1CH/O07oBsIRfRvzNBrcUjAoqVeE7wCSjx5AvmahUBZbyVr7hyr2jyBsVPDufq3lsmD49ErkKV0b9MkGtqzFMpLU6q4fLORmGGRBkn5564U96hF9G+Ym44mxzI34VUxkdRqWKtSbtlf06c0rSdofOpbyBqgAmh3RctXcEAPDITG2C9/7nrKGDF06lW7432XdHiror/ArG01pNxF8yq4jHCDGfGb6dwMLPMExkkRF/fXLX8vjtRVu2oO8cTngifktM03bED7jRuB9uL34/q8cSfmk7TWV0XDCVQUZXGxZyffXhWVy7fxR7x1It35v062Xi1mnLELeEvz7i92vQ1iks/AzD1GCYVdzz4+5OvatUBf752NmWc3LrmVsrQSE0jDFMeco5pdjvHErUdOkErDsDWRffrD1zswEmchiLrOGfGtIRUwiH9wzXVPY8dW4NT53L4nVX7wn13qRfLyN+afloMUv4l+uqepo1aGsXFn6GYWq48+EzeOun7vedCdsp33tmHu/63IN48NRKW8+bz5YwkdEbLI6UTznn1LDeYPWkddW5aDRL8Pq1ZJZkEtb4xYW6fMN1B8dx7MwqHp+17J6vPjyLmEJ4zZW7Q703J+I3a4Vfj8carR6febsbgYWfYZga5OzaMKtdwzKzbK1yXV0PNxRF4p216yVV5/En4zEMJ9xpWbIzZzIew+RQCKtHevw+Vk+6LuKXdxDveMkhjKc1vP+OYzArVdz58Cxe8rzJpiMXvciIv1Rv9agKJtIaciXT2Vbymbe7EVj4GYap4fh8HkD7It2Ms6uW8GfbnF87ly35Cn9SU7FerqBaFciVTGQSqrWat87qSesqxlMaiID5JlaPM33Lx+oZ8nj8KS3mdOwcTWn4wM2X4ZGZVfyvLz2CMyvreN3V06HfW31y17F6VAVjdssHGfUXy5Wu1fADLPwMw9QhSxRXCuWu7fPsilWd4+2lE4a5uj49EmcYS7mCXKmCjK4iravIG9bFIO+p41djCsZSWovkrvT4/SP+klnFudViQ3XRa6+axssumsKdD89CVxX8zOW7Qr836dkXG8o5FafXjyP8ZqVrq3YBFn6Gwad/cAL3n1jq92FEgpJZcYaDdFX4Vy3hl83OwlCpCizmSg1iC1jVOoDl5eeKZaT1mDuLt1zBulEBkdviYCKtNff4i5ZdpPmIq+zQeWIh32DjEBH++PWHkYgruPGync5jw6AHRPy6GsO4nZeQwl8qVzniZ5husZw38IdffwKf+9HJfh9KRzx0ahnv/vyDDStIO+XUYgFyV/W94jeCtHpybUT8KwW7bt5nFWxSk8NYKsh7In7AGsaSN0ykNdVp5TyZ0Zu2bVhbL/tW9ACu8D+3mK9pFifZN57C13/jpfij1x8O/d4AbzmnLfxl1+Mf94n4WfgZpkvc/eM5VEW4QR2Pz67ipr/4HrLF7kXCG+V7zyzgG4+edWrMN8pxz0rUlfXuvE8hhBvxt+Hxr5fdRVj1uOMXTcvj19WaoSnrRsXp2w9YLRQW883LOf0Su4A7hatgVHzvPgDgeTsyoWfySuqTu0aldgEX4LaAKJarbPUwTLe460lrhN5CtnV0+9iZVTxxdg2ztl8dBeRFyNu7fSPIxO5kRu+a1bNcKDs2hrevTivc8sZGmfIOY5HCLy8Q+VIFBaPi2EFAbflnPQu5Eh44ueybRAbcnvxA4wrijdBQzll2rZ7RZBwK9S65y905mYGlXKniu/ZCpTARv1xhuV5uL0HZS2Q1Sv1w7k45Pp/DruEEdg7rXYv4Zz0DS3LtCH+TVsTpGqvHRFpXHY8/b5goGKZjB1n7iDmlkV7KlSre9bkHsVIo4/03Xep7HBndFdywpZphkELulnO6VT2KQhhLaViy7TYu52SYLvEfzy0hWzRx6e5hLBUMmJVq08dLwS8Y7ZUkSmRk2U1k/Xm3hP/Z+TwumEpjJKVhtUsev7R54jFqT/jN4MZkKSfiN5H1lHMC1l1FfcSvqYpTNePlD772BO4/sYQ/ueVKHN4z4nscGd21gLoa8au1K3cNTx0/YF1kTtuJ9iIv4GKY7nDXk3PQVAWvv3oaQrQWT2kV1A/PCMvffv8E3vLJ+1DtUiIWcCP+Zu0IwiKEwPH5HC6cymA0Ge9axH/OTuwemky3KfzBEb+0elbXyzDMKjKam9zNlUzk6zx+XW3sff/D4wv4+/tO4p0vu6Bpm4W0J+LvpvCrMQUxhXxW7lrv9+UXT+He44tYzJVQ4pYNDNMd/v2pObzwggkcmLAaarVaqSr/QNeN5ncGQSzlDBTLVax2SVABOB0nl7rg8c/nSsgWTVwwlcZYKt41j392tQhVIRyYSHfN45cRv/ydZRJucrdgVLBuV/VINFWBWRU11U8ymn7rCw80PY6hHkX8gBX1F+utnpj1ft9w7R6YVYE7H56FUalykzaG2SjH53M4sZDHT1+6I1QTLwAoGhuzerJ1A7u7QXa9e1bPs3Zi98KpDEZSGtaK5a6UiZ5dWcfO4QSGEmpbK3dLTXrQy0Tu3Jp1LtO6ipT0+Esm8qWKc3Hw7sNr95Sc4SbNBdUb8YcZot4O3vGLhllFTCGotvBfsmsYl+0exj/+xykA/hfATgncExEliGjKZ/sOImpcSscwW4jvPGVV87zyEo/wtyiJlB5/p1aPd4pTt5ARf5jRgq2QpZwXTKUxmoxDCHSldPXsahHTowkM6aozwjAMzaweJ+K3z+WQ7nr8uZKJ9XKt1SMXZtUIf8g5tnL84nBC7WplDSCF3+3VU38sb7h2D3583vq9bFbE/3EAL/XZfiOAP+/aETBMD1grlpu2Fp5ZXsdwQsXesRQmh1p3bwQ2XtWz5vjxXYz4nXLO7kT8ibiC6ZEkRlOWvdENu+fsahG7RpJOs7OwrZm9K1nriccUxGPkCH9aVxFTCMm41bytYJiO52/tQ067cn93bt18a0HN6GrXbR7AiuK9Hn/9yuHXXj0N2Zh0sxZwvUQIcUf9RiHE5wC8rGtHwDA94HP3ncLbPn0/VgOEK1s0nb4saS2GZDzWMhJ3q3o6E/6cLdLdivhLZsURx25YPcfnczg0mYGikCv8G8xHCCFwbrWI6ZEEMgkVZlU0JFmD8A4m8SMZj2HOntAla+3TegzZYhnFchXJeGPEX6qJ+K22DvFY66lWaV3taimnJKHGHEvLMBsXae0YSuBlF1nGy2aVczY7G5wbYHrCp75/Aj88vtDWcx46tYz7nl2s2XZiIQchgv30XMldok9EmBxq3ssFcC2eTiN+P6vHMKv47L3PtSwlbbY/VaGmDcjC8ux83hkZOJK0vOyNtm1YzBswKlXsHknUrKwNQysrJq2rzrmU+7a2GfbXtVU9QJ3w20Ir2zo04xUXTeEVF+8IddztkIi71UbW8TRG9W+4di8A1FzINkozAZ8jouvrNxLRdQC6O56HYWw+9m9P47M/bK9vzp9/+xl84CuP1WyTjcaCBDFbNGsaak1m9NbJXenxdxjx+wn/d388jw/e+Tgeqhvh187+9o2nsFY0fevUw2KYVZxeLuCCSUv4ZcS/0Qok2ZVz10jSOd9hK3uaWT2AVdIpu3A6wq+pzsW+dgGXj8cfILR+/P7rDuPXX3FhqMe2gze5WzIrvk3ifvbwLnzo5y7DSy9qSLl2TLOVu78F4ItE9BkAD9jbjgB4K4Bbu3YEDGOTt+uvn13ItX6wh6JRwcnFAipV4UxqOr1k1Y4HWSC5kun0QwFqF8sEsRGrp1ypOs/33oXMLFuvGdROoBmylfDBiRROLOSxXDCwc7izuovZlXUIYV1EAGA02R2PXzZnmx5NOAu5wlb2SKvHTwwB1FTtpHXX6jljD31J+1T1eD3+IKHdTHRVcc6Hn9UDWPmM//riQ1193cB3LYS4H8D1sCyfX7b/EYAbhBA/6upRMAzcSPg5W8TDUqpUYVSqTmsAw6xi1hacoKRnzifib+W9byS5621O5n0dOZkqrO/tRQrGQTtK30hlz2n7AiSHhI90Tfgtsd89ksRQBxF/PEYNYxclqbj7+5Min9bdiD8Voqqnm43POqE24t+842naq0cIMQfgQ5tyJMzAI/9gDdMScRl9tkImx55dyGPfeApn7OgVCI74syWzpg3v1JDutG2QddT1bGTlrhRpXVVqcglnHOHvZJ8y4reEfyMJXnkB2jeeBGCVMA7pKlY2OIVrdnUd8RhhIq1htgOPv5kVI+v2k/GY8ztLayrKFeuXn/Kxevw8/n6SiMdqmrRt1h0IJ2mZyCAX4wCWiIdFluWdsOvQT3ksm0Crpy7in8poVtuGJsnM4gasHtlT59BkGot5ty/QGfsuRSYyO9mnE/FvYPXu6aUCYgphl8cqGknFA6uiwnJutYhdIwkoCjmVN6GFv8XUqZQnypd4E7qtF3B1t/9NJyTiSl0d/+YcDws/ExlkaR7gingY5B/zCfticWrR+j+jq76VOqbtt3ubb7mLuFoLfydWj4z4L5hK1/QFkh7/RqyeQ12K+KdHEzV3O6OpjffrObtSxO5h6y6i7aqeFhF50rZ6vHdu3ouA3wKuhoi/iyWSnaCr/bF62noVIlKIaLhXB8MMNvPZElSFkNFVR8TDIP+YTyxaInpqqQBdVXDRzoyvGMq5r94+63IRV1D5pxDCEfwwidi1YrnGy5Zid8FkBoA1SzZfMrFsR9SdWD1r62UQWYlThTbm8c8sF7B3tNZaG01qGy7nPLu2jt2j1l1EJ1U9epMSRhnde6N8b3+etK/V403uRsPqkZ9fw2cBV69o+SpE9HkiGiaiNIAnADxNRL/V+0NjBo25bAmTGR0XTKXbs3qciN+1evaNpzCR0X2FX/bMGapL7gLBbRuMStUZSRjG4/+1v38AH7zzcfc1PVYPYK3ePePpU99JxL9m21XOMPEWEf+jMys4t+o/ROb08rrj70tGuhDxL2QNZ1xhSouBKPwUrlK5udUjI/qM7h/xt0zutlHO2SsScatddNVe2BYlq+cyIcQagNcD+GcA+wH80kZelIj+BxE9RkSPE9F7NrIvZvswly1hx7COCybTTsOwMMgobmZ53R4Wvo794ylMpP3FUEbffpOVghZxFT0efBiP/9xa0bkQAR5bxl4gNZ8tOYldAB3V4GeLJobt1ccTGa1lh853fvYBfPzfn2nYXixXMJ8tORU9ktGk6/EbZhX/64uP4Klza6GPr1iuYL1cwZhdNktEyGhq6Lm7rSJyWdXjFX7v0JSUb8uG2pW7/Y743TLT6qZaT2FeJU5EcVjCf6cQogyg45Z9RHQYwK/AKhW9CsDNRPRTne6P2T7MZ0vYMaTj0GQGs6vroSJrIQQMs4p940kIAZxcLOD0UsES/oyGpbzR0P9eRpw1kaIWQyKuNBF+t6Y8jMdfKldrLjoy4r/QtnrmcyXH3wc6tHqK7urj8bTW0uNfLhiYW2uM+OsreiTS4xdC4NiZVXz5wRl89+nwazeXbZtozDOLNq2ryJXC3UW0SnamfCJ+byWPX8sG7wV2M62VIGQbhmK5Yq0rCKgo6zZhXuWvATwHIA3gHiI6ACD8Zb+RSwHcJ4QoCCFMAN8F8PMb2B+zTZjPFjE1pOOQnQB9brF11G9WBaoCuHinlXp68OQyciUT+8dTGE/rqFSFU/0iyfpE/ETUtJZf+voTaS3UBalkVrCU8wq/CU1VMJKKI2O3GphZWYcWU5CMxzqq6skWy27En9abevxmxYoo/e6AZupq+CWjSQ2VqkCuZOJBe3JYOwnk5bx13sdSbhI9k1CdHEsrWkXAKd2vqsf6OhFXaur//RdwRcHqcefuRiriF0J8XAixRwhxk7A4CeCVG3jNxwC8jIgmiCgF4CYA++ofRETvJKKjRHR0fp47RGx3zIolSlNDCadtwIkQdo+M4C7ZNQQAzgxdafUAjYu4ZMTv9fiB5m0bZK31WEpDuSJQbtFbp1iuIlsyHaFZK5oYti80kxnNsXqmRxNIarGOq3q8EX8zj1+2NvAT7tMy4q8T/hFPh84HT3Ug/DLiT9dG/Nm26vhbl3NmEo3lnN7IH3AbsTWUc/a5qkdG/OtGxV65GxGPn4h2EtHfEtE37e8vA/C2Tl9QCPEkgP8D4FsA/gXAIwAaPglCiNuEEEeEEEemprrXo4KJJot5A0LAtnos4Q+T4JWCOZnRMJnR8P2fWA3e9k+knJYM9WLl5/Fb+9Adq+fup+fwD/e5PYNkxC/32crukXcF8rWzxbLTDXRqyHqdmeV17B1L2WMBO7N6hpOux7+6Xg68IMlKmiWfC9vMcgFaTMGOurbDsm3DcsFwhH+5jSofP6tnSFfbqOppbvXIcs6M1hjxexO7gHVHVz9+MRIrd+33513gtxmEeZXPAPhXANP29z8G8J6NvKgQ4m+FENcKIV4GYAlAY8aJGSjk4q2pIR1pXcXOYT1UgldGcHo8hkOTabdx2Zgr/PUWiJ/HL197IVfCudUifuPzD+Hjd7kfSyn0snlZs0ZtZqUK084ryNfOeVYKTw1ZltKZlXXsGU36zoMNgzfil3c3QcIsp4Z570IkM0vr2DOWhFLXGmHUFuwnz67hvP37ac/qkRG/a/Wk9Vj4qp5WyV2fiD8TIPyA5fM3rtyNhtUj7cgoCf+kEOKLAKoAYPvynbUntCGiHfb/+wG8AcAXNrI/Zuszn7OSjjLqvGAyU1MVE4TTyCumOHcKO4Z0JLWYU6JZv6JVWg3pOjtgKmPZJb97x6PIlsyahUbSg5cXk2aVPV5xkXcQXpGeyuiYXSliPlvCnrEkdLV9j18IUWf12O81wKryVtJI710ys1zA3rFk/VOci9y/29PKDk2mnXUHYZCPHU26EX9Gj7e3gKuJFZP28fil4CfrfreA5fPL340QAkal/xG/fH9r69GL+PNENAG7koeIXgBgdYOv+2UiegLA1wC8WwixvMH9MVscGfHvsFsGHApZy+9G/AoO2RUz++0ePzLSrLc3ZLuG+gh3ckiHEMB3np7HntEkCkbFaRYnI/4wVo9X+GusHt21euTz944locfbt3rkscnkbpCt5TzeI7b1F8LTtuVUj7R6vv/MAhJxBS+8cKKtiH8pb2BIV2sqZzJ6rI1ePc2tHnmxm/IMSHHbMzc+T1cV5/PSbJD7ZiLfn2x/vVlVRk2btNn8TwBfBXAhEf0AwBSAWzbyokIIv5GOzAAzZ1fTTGYsAbtgMo2VQhnLeaMmOViP/AO2In5LvKTw62oMQ7ramNwtlRtsHuu1LQG5Zv8oXnPFbvzRN55ErmhiJBV3PHvpVzcTfm/Vj4zAs0XTsSS8k5w6tXqkpSXzBvK8BXYj9YitV7zzJRNLecM34pf5g7xRwfWHxrFjSMfqerlpIzsvKwUDox6bB7Ci83zJGr/YagBKK6vn0GQaX3n3i3HlnhFnm0zq+lk93lxKq17/m4VM7rpWT0SSu0KIBwG8HMCLAPwqgMuFEI/2+sCYwWI+W8JoKu588GXHyVYlnaUaj9+K+L1dPcczjfXtuZLZkNgFgKv3jeK6g2P401uuckRPrvKVQi8vQs3aNniFfyHvY/V4kqh7x1M1FkRY5LoAb1UPACwFrEPwDjn3ng+3hr8x4k/EY04t/PMPjDmvEXY171KhjPFU7UU77PhFIURL4Qes35n3zk1TFWgxpaGqR/7MjfjtsY79Tu7GayP+yFg9dsnl+wC8RwjxGICDRHRzz4+M2bb4DdueyxZrqkqkMLdaJev1+C+YSuN1V0/jZy7f6fzcKnOs8/jrOnNKpkeT+NKvvQjP25FxSj1lZO2t4/d+739MHqsnZzi18N6qHgCIKYSdQ3pHVT0yQpQXqNGUBqLgiN9bO+/NA7g1/I0Rv7Vfa//X7h9z7nbC2j0rBcNJEEvCNmpzBqF3MG5wOBmvadwm8d5ZyZxK/xdw9cfqCfMqnwZgAHih/f0MgD/q2REx25pf+/sH8OGvPt6wfS5bqomEnZWWLerlvR5/PKbgL269BpdPu7f+fgubcnW9+P2QIi0FSgqGFMLQVk/ecPYxXBfx7x6xumHqcaXt5O6aY/VY+4wp1LRfTz7A6pkJqOGXyIEs1+wfbZlHqGcpb9RMOQM8wt+isse1YtoXwr/8hWvway9vHJOo+Xn8/RZ+VSZ3I2b1ALhQCPEnAMoAIIRYR/NB7AwTiKIA33zsXEMbBatdg9sLXi5db9XDxvB4/H5M+LQyqO/F74e845CWyrpRgUJwkqlhIn5NVbCYNxpsmQk7Kbln1IqyO7N6ai8m1n61wBnDec8di/ficNruZCpzBPWMpzUcmEhhMqM7Ef9y6Ii/7FwoJemQEX+rQevNeMEFE77WlXWe3fm2cls/abB6NinZHCa5axBREm5Vz4UAOp/4wAw0r75kJ/752Dk8NruKK/eOArCsnzm7T49ERvytVsi2itzGMxqWC0ZNMjFXCiH89VZPuYJkPOYkDcNE/NMjCSzmSo7IybsIzRZaNwndgdVjC4W8EAFW0jionDNfMpHWYg3N3GQn06BE62//p0uci+uEfXFoNqxGYphVa65xndUzFFb4eyDMmqpgZT1qVT0yuWudj83q1RNG+D8Ea4XtPiL6HIAXw5q/yzBt88pLdkAh4NtPnHeEf61owjCrtVZPmxF/kEBMpK0WC2tF07EtckX/5K6X4USt8BfLFSTiMSdpGKacc3o0iYdOrXgqcNzX/OtfOoJdI9YdTv3CojDUV/UAloX0yMyK7+PzJRMpXbVyHp6Lw6mlAg40GXF59b5R52sZvYeJ+GUf/9E6q0dG/K1W7/ZCmHXVtdQ2ckfRTdSYAlUhZDc54m/6KkSkABiDtcjql2EttDoihLi750fGbEvG0xqef2AM335yztk2b0/e8gp/3Keboh9eWyXo9QA4Fki1KpAzzIY+PfXUjwlct4VfCkWYqp7p0STWyxVnjYJXpJ9/YMxj9Shtt2VeK5ahKuSUAwLNB8bnjQoyuoqJtDujQAiBU0sF7J8IN9tYV2PI6CqW8q2reuTiLb+qHqC3Vk8Quqo4OSMnedxnqwew7J7VKHn8QogqgP8mhFgUQnxDCPF1IcTCphwZs2356Ut34omza5i1B5E4i7d8PP5WVo/Roixvwq6Zl2JXKFcgRGOfnnqS8RhiCjn+fLFcQVKLQbHFtmnEb4uWFHZZkhqUUO60nHM4Ga+xaKaGdBSMim80nS+ZSOuxmmZu87kSCkalacRfz1g6HqpfjzzfY3Uef9iqnl5ZPW7EH41yTsCq5Y9iVc+3iOg3iWgfEY3Lfz0/Mmbb8upLrXLLu548D8Add+hX1dNKEFtF/PUdOt0+PXHfx0uIrBGQOcfqqTrRdUpTWyR3rZ85wm+vQA66y9BVBZWqcAawh8G7LkDSbJhMvmQipVlWj2zmdsoeVXnAXjMRhvFU677/gGv11C++24yqniB0NeZE+nL/iT57/IB1XLK3U2Tq+AG8HcC7AdwD4AH739FeHhSzvblwKo2DEynH7jlvDwfZMdzo8ZcrzWf+GK2Su3UliHIISKuIH7AidNnXZ92oOIuZkvFYi+Su6/ED3ojf/2Ijfd12ov619XKg8PvZPXnDSmjL6p3lgoGTtvCHtXoAS8jDCP+ST2dOwB2/2A+P34r4a1fuarEoWD3ue9ws4W/56RdCHNqMA2EGByLCT1+6E5+99yTe9+VHccdDZ7BrOFETEftNTPKjZFahEAJbCNR7/E5StIXHD1jRqbeqRwptIq6Eivin7SHjzy0WGvx4L97xe2nd9yENeMcuSqSo+wl/oVRBaiLm9LdZyhs4uVQAUfDiLT/G0xqeOd+6ed6KbNBWZ/XI8YutevK7Vkz3hNnr8TtWUgQi/oTPpLBe0/LTT0Rv8Nm8CuCYEGLO52cM05IbL9uJT37/BP7poTO45fl78asvu6DGr44pBIUAo9K8zNHqsBgsDom4lZB0rJ6AXvx+DCfiHqun4kTUSa11xK8QsNNuODefLWEsFQ8smXTnwYYv6cwWTRycrI3Um1k9soTVbe1g4NRiHtMjybbEdTylhfb4U1qsRtQksl9PM3ph9WiqgnJFWIPNI1LVA9QJf4TKOd8Ba9Xud+zvXwHgPgAXEdEfCCH+vkfHxmxjrj80js/81+tweM9ITdMyL/IPtRmlcqVllOSdRxvUi9+PTELFnF1xVCy7Vk8q3tzjL9pdJdO6ikRcQbFcDbR5AI/V08bq3TXPYBfJRFqHQgFWT8lEWledWvxFO+I/0IbNA1hWT8GoOOWtQSwXjAabR5JJqCGSu92vupH7MuwxlN3ef6fIi4+uKi0b13WLMJeXKoBLhRBvFEK8EcBlsBZw3QDgd3p5cMz2hYjwiot3BIo+YEU/Lev4Q/RU3z2SwOkly8925u2GEP6hRG1yVwp/Qouh0KKOX9o6cpVusxYRXqsnLH5WT0whjKd1J1kuqVYFCuUK0lqsJudxarF94ZfPbxX1W11V/S921sD1cD2Yuu3xA9YF1unxFKGIfzOPJcwrHRRCnPd8PwfgIiHEEuw2DgzTCzSPJxtEqVxt+Qdz8a4h/Ph8DkIId95uCKun3uOXYp6MK00ncHmjYem7Nxf+9qwet+mbX2tpDfPZWlFet0tY07qKMbuZ26mlAhbzBvaPh6/oARC6UdtyoRwY8Q/pKnLF5tLRqzp+AChVrMHm8RjVDGTvF/JztZl3H2Gsnu8R0dcBfMn+/hYA9xBRGsBKrw6MYcJE/KUQEf/Fu4aQK5mYWV53LIZ0qIg/7lb1lCtIaO4g71Yrd+UxyQi5qdXTZsTf7OI1NdQY8cuWzCldRUwhjCbjePj0CgB0HvG3WMS1XDCclhT1pPVY4EIzSS+sGG/Ev5mDzVsh5+5uZr4hjPC/G9bK3ZfAas72dwC+LKzeuq/s4bExA048xIpWK+Jv/gd8ya4hAMDT57LIlUwk7E6erRhKqDDMKorlCgzTY/XEY03bRXsjfrmArFkVUbse/8q63Q7BJ6KeGmqcVSxbMmfsUYXjaQ3HzlhD9ILEOQinSirfXLiXfTpzSoYScaetdBC9sGKksFoefyUSiV3AbT0dKeEXQggiOgpgVQjxbbs/fwZAtudHxww0WkxpvXK30trquWinLfzns3Yv/uaLtyQyopbRacJTx18MGfFPpLtv9TjtEHw8dDnI3duUTlbQyBnDE2kdx+2LQ+cRf7DVY1aqWCuaDaWc3n0s5Y2mU7h6YcXoXo+/3P95uxJp9UTK4yeiXwFwO4C/tjftAfCVHh4TwwCo7Z8ehDWXtfnHeCgRx96xJJ6yI/4w/j7gJoCldeJU9djlnH4DZQC7qseJ+Ltv9UjR9Y34MzoMW3gl+Tp7S4r3eFprelx+jCTjILKmawUhJ3QFRfzjaQ0ls9p8YH25+1ZMfVVPFBK7gBtQdDJ0plPCvPN3w+rIuQYAQohnAOzo5UExDADEY62Tu2GqegDL7nnq7BpyRf95u35IUZQRv7NyV4uhUhWBx1Y0q67VE6KqR2s74reEv74BGuC/eld6/I7w2xejdqN9AE6OoFnE3+zCBLQeDA+gJ1aM6/FX7P1Hw+N3yjk3qYYfCCf8JSGE8xsiIhV2b36G6SVhIn4jxFxWwErwPruQx1LeCC388nFyELz04qWoFw3/Y/PehYyHivjb8/jdBmj+ET9QJ/y2x5/W5MXIFv42/X3JWFpr2pM/qDOnpL5/kh9h5u22S63HX43Eql3AG/FHS/i/S0TvB5AkohthVfd8rbeHxTC1S+yDCHvLfvGuYVSqAk+ezYZatQs0evzeXj1AcE/+kifil0Ic5HcDXo8/ZHK3UEZMocCqHgA1lT1BVs/+NpqzeRlPac0jftmLv4nHDzTPE1jC3N2IXIuqx+9ZwLVZhHml9wGYB3AMwK8C+GcAv9fLg2IYwLJ6WrdlDucFX2pX9hiVaqg+PYCP8Guuxw80EX5PxH/59DD+7D9fhVddEuyOSoEL25N/qWBgLBWH4pP4lAviFmqsHhnx1wp/pxH/eItGbVLQm3n8QIuIP0Tupl1qPf7oWD39WMAVpqqnSkRfAfAVIcR87w+JYSxC1fGblVD9TQ5Opq39VaqhI34nuVtX1SP/Lxj+bQeKnpW7RIQ3XLu36eu0XdWTNwL985FkHPEYBUT81nFfOJUBEXD5nuFQr1fPeFoLnPQFuFZP0AIu1+MPLgnthdXjzaWUzCrGUhGJ+J1yzggkd8niw0S0AOApAE8T0TwRfXDTjo4ZaOIhevUYIb3aeEzBhTsyAMK1awDcRm71VT0y8g8q6SyWK86inDCodkO60FU9BSPQP1cUapjElS+Z0FXF6WB6eM8IHvy9G3HJrs6Ef8xTjunHUr6ERFxxzlM9GV2FFlMaIn5v/6NeROSOx29WQ39uNgN35W40rJ73wKrmuU4IMSGEGIfVn+fFRPTezTg4ZrAJF/FXQ3c0lAu5wkb8uhqDpiqObZKI11k9QcndNkWFiNqawrWcLzfNGUwN6TUdOvOG2bBSuX5ASjuMp6w5xkGtlc+uFrFrOOH7M8B6v+NpDUue2b/3/HgeV//BvznH3Yvkq3e4TylCK3f1iPXqeSuANwshTsgNQohnAfyi/TOG6SlhhpC3E7ldbAt/WI8fsIauu1aP7NUTbPWUK1VUqqKtiB+wKjpKTRaFeVkuBK+KBRpn7+ZLFcfm6QaTQ9ZrLwS0XTi7WsTukeY9/sfTte2dH5tdRcms4pTdTK8XydeaBVwhLcLNQI9YcjfuN1/X9vnbW/XBMB2gxahpcrdSFTCrIvQUpXYjfsCyJWRlUbLO4/dL7nY6OUoPcZEDrAHpy4Vgjx+wKonqrR6Z2O0GTgI555+cPbdaxO7R4IgfQM3sXwA4u2K1v16099kbq6d2AVd0rJ4IefwAmrXfaz2JgWE2SKs6fqNNkT1ycByvuWI3rjsYfmS0t/5eetbNPH65rVmvej/CWj15o4JyRfi2a5BMDelYzBuo2nNc/ayejeAKf2PEX6kKnFsrYvdIa+H3VgbNrqwDcCel9SK5G49ZVVClciVi5ZzRquq5iojWfLYTgOa/VYbpAq3KOZ1GXiFv2TO6iv/7lmvbOgZvIlj+gaYcq6dJxN/mH7EV8be2elqtigUs4a9UrTuDiYyOXKmCkWT3btKbTfqaz5ZQqYpQVo/X459dLdbssxcRuZVLUWyPP0rlnJtv9QQKvxAiGmeFGVg0VYFZtUbl+dWstxvxd4Ks5ddUxTmGZJM6/o4j/rgSauVus3YNEhmRz+dKmMjoKJRMTLeIwNthLKUFTvqaXbUi9+kWVs9EWkO2ZDoCfNZ+nrSPSuXeCLOmKsgbJqoiGmMXAa/VEw2Pn2H6ijNwPSDql9F1L5N0Mh+Q9Ai5NSIPvsNYpPC3K1pBVs9iroQ//9aPYdrnwGnX0MLqAVxhlmMXu4Wc9OUX8Z+zI/ddwy0ifruVxUqhjHWj4gxnl75/L6wewDrPcrhOVDz+HUM6XnvVNF5w4cSmvWY03jnD+CAFvaXw9zBSkuMNEx6RICIkA3ryd5rc1WL+Vs+3njiPv7jrGTw+a7muKy0WRwHA3jFLdE8sWK2X80bF6dPTLfwmfQGuV98q4pd3LIs5w7lLAKxKISFED4VfcYU/IlaPGlPw8Tdf0/G6ik5g4Wcii+ZZcOOHM5e1h3/A0uNP1lk3yXisudXTSTmnz/uULY6fW7REfKlFOwTAmjE8mdHwyOlVCCG6HvEDjWsFJGdXi0jGYy1zCt4OnbKiZzQVx2K+5Fzoe9Gm2BL+svP1oDK475yJPHJKVlCC1+gwkdoO0uOv9+wTAcK/oXJOH49fevoyel8uGFAIDYPWvRARrtw7ikdnVmBUqjCrovvCX1cyKjm3alX0BA1Ykcg5BYv5khPxX7FnBIs5o+MEeRg0VcHauul8PagM7jtnIo9j9QRG/L0X/kyA8Ke0WE2LAeeYOo341Ziv1bNqWzvPeYR/NKX5Jru9XLl3BD+Zz2FuzRLnrls9dsRf37ZhdnW9ZQ0/AIzbcwqW8oZjDx3eM4KlguGc195ZPTLij4bV0w9Y+JnI0srqMTbB45d1/A1WjxZk9VjHlOjSAi4n4l+0VrS2atcguXLvCIQA7j+xBCDccPl2mMxYU7RydW0bzq60XrULAKPJOBRyrZ7JjI7pkQSEsOwioDfCrKmKM52MrR6GiSDxFsld1+rpXeQm2zvUNxwLGrju5B06KOf0u8DJZO7JRTfib1bKKbly7ygA4IfHFwF0X/jdWn43wWtWqpjLtl68BVjN5MZS1iKu2dV1TI8mnMH08g6gF1U3uhpzLlZRqerpB31550T0XiJ6nIgeI6IvEBEvCGMa0Fsmd3sf8btWT+1rXDiVxtHnlvDFo6drtjsRf9sLuPzLOaXwrxTKWCkYWGrSktnLZEbHntEk7nu2N8I/6TPpay5bQlUgVMQPuF0+z9p5ATmZyxH+Hnn8ErZ6NhEi2gPgvwM4IoQ4DCAG4NbNPg4m+rjJXf/2v0ald16wJCi5+8GbL8eLnzeJ3779UXzmB04fw84j/oCVuyvrhiOIJxbyWCmUm7Zr8HLFnhGcsUW0++Wcjat35SKsMBE/4PbrObuyjt0jSUzadxEzy1L4e1PV4/f1oNGvd67CGuWoAkgBmO3TcTARpmU5Z3kTIv6gck4thk++7Qj+n8t34sNfewLfeWoOwMYi/nJFoFJ1L3JWQ7Yyrtk/CsAq6bSmb4VrqXzlvhHn615F/LXCb3nzYZK7gLV69+RiHnmjgj2jSUym66yeXkf8bPVsHkKIMwA+CuAUgLMAVoUQ/7bZx8FEH3flrn8PG6feexOSu34tGHQ1hr+49RoAwLEzqwCsiD+mkDP0JCxShLwXuWLZGhhyxZ5RKAQ8MbsGw6yG7qV/le3zA+GHz4RlPG21bfC2Zpb1+GGtnvG0hvN21dHu0QSGkyriMXLKO3vj8bPVA/TH6hkD8DoAhwBMA0gT0S/6PO6dRHSUiI7Oz/PEx0FEdlM0TH+rZ7MiflWhQOFMxGMYTcWdyLdYrrYd7QP+4xdlRc+OYR3To0k8dGoFADAWoqoHsMojJakuWz2ybYN3xOPs6jpSWgzDIdteT3guYLtHkiAiTKR1zK70tqpHwlbP5vLTAE4IIeaFEGUAdwB4Uf2DhBC3CSGOCCGOTE1NbfpBMv3HSe4GVfVUel/VE1MIf/vL1+EtL9gf+BjvYqZiudJ2gzbAfQ/eBK/bniGOQ5Np564irNUzkrSeB3Tf6gEa2zaEXbwl8d65yBYPExm3XXOvevVIeAHX5nIKwAuIKEXWJ+TVAJ7sw3EwEcdJ7gZ6/BX7ceGEplNeftEUdgwF+9ZTQ67wd9pjxjsdSrJiR/wjSQ0HJ9LORaFZu4Z6rtw7gphCPRHR+rYNs6tFTI+Gs3kA933EFHLOryzpBDji7yXdDwNaIIT4ERHdDuBBACaAhwDcttnHwUSflt05K1VoqhI6wuwVU0M6Hj69AmADEX+80eqRfXrG0nEctCN3oHkv/nre8ZJDODw90pNzNJXRnVYSAHB2ZR0XXRT+7nzCTubuGk4gZq9EnvRc1Njj7x2bLvwAIIT4EIAP9eO1ma2D1qJXT1SmKHln3JbMakcWgnyvXqtHevyjSQ0HJ1LO9vYi/lFnMVe3mbTvdIQQKFcE5nMl7O4g4veWf8qSTqC3Vg9R7+8Uo0xfhJ9hwhBv1bKhEg3hnxrSUTAqyJfMDUT80uP3RPy2xz+aciN+InR1mtZG8LZtWCmUIQTaGvgiG7V5LxbehG8vrR49AneK/YSFn4ksflGwFyvi7//t+pSnpr1UrrbdpwcI9vgTcQWJeAz7xlJWV85k3LFF+o136Po9P7Yq77yVRK2Qk7z2eIXf4/H3IvmqO8Lf/89NP2HhZyJLK6vHqHRmq3Qb78Srkllpy4OXuOWctVU9o0lrX5qqYO9YKjKiD9S+73+47ySu3DvSlvBrqoJP/OLza6yoSfsuIB6jnrxXb8Q/yLDwM5FFUQiqQk1W7lYi8Qfs7VtT7Djib7R6lgu1nTiv2jeKQl03zH4i3/c/HzuLZ+Zy+JM3Xtn2Pn7m8l2+++xVRO5E/AO8ahdg4WcijqYqWyfiz5Wc4eHt4lb1uO91dd2oEf4/vaV9Ye0lUqQ/f/8pDCdU/NxV0xvep/T9e3VBl/vt5ZzmrcBgv3sm8sRj/u2KASvpG4WI39u+oPOI36+qp1yzWCsRj3WUOO4V8n0bZhVvfP7ehtbVnSBLPHsn/LGa/weV/v/VMEwTNFVpOmw9ChF/TCFMZKz2BUWzuyt3wwxd6ReybQMAvOWGA13Zp6YqGE6oPZm3K/cPsNXDVg8TabSYEtirxzCrESpttGraO11b4Fg99mpkIYRt9bSfKN5MDk6kcOnuITxvR6Zr+5zM6D27oOuc3AXAws9EnOYRfyUyXq1s29B5xF9r9eSNCsoVgdGIXNiC+Ju3HoHa5YVQExkt0N7bKBqXcwJg4WcijhZTAnv1GGY1MrfsUxkdT55dgxCdRZP1axZkn56wDdn6RdgW0e3wG6/6qZ4Jv+vxR+Nz0y9Y+JlIE1epuccfsYgf8O/d3woiqpnCJVftjkTY4+8VL2uj30+7uB7/YEf80firYZgAtBZVPVFI7gLuwiOg82hSVxVn5a7bkjnaEf9Wgz1+i8F+90zkiceaV/VExaud8jYX6zCa1DwD150GbQMY8fcSGShEJWDoF4P97pnIo6lbI+L3Cn+ntfY1Vs+626CN6R4c8VsM9rtnIo8esHK3WhWR6c4JADu60E5YjytucjfvtmRmugdX9VhE46+GYQIIWrkr7Z/IRPwZtx1x5xF/zPX418tIa7HIvL/tghZTsHNYx77x8HMDtiNc1cNEmqA6fhkZRyXiH06qViJ6A3chXqtnuRD9xVtbESLCPb/9SsSVaHxu+sVgv3sm8sQD6viNiAk/ETmVPRvz+K33tRrxdg1bGV2NQYlQe+t+EI2/GoYJIDjityLjKHm1MsHbSZM2wKoG8lb1sPAzvYKFn4k0QXX8cluUPHAp/J1ejKw6freqhxO7TK+Izl8Nw/iwVTx+wO1P33HErypYyhu468nzWMxxxM/0juj81TCMD4MU8R+aTGMuW8I7/u4oVtfLNWsDGKabcFUPE2niMQVVAVSqomYGq7wLiJLH/6ILJ/Hw6RUMJzr7s/qfN16EX3zBAcyurGM+W8INF0x0+QgZxoKFn4k0MqI3zGrNhCdZ7x6liP+FF07ghRd2LtZEhJ3DCewcTrR+MMNsABZ+JtI4wl+pIokYbn9gBpMZDZWqNZwlSh4/w2wVWPiZSKPZQz4MswrDrOI3v/QIAGB6xIqKoxTxM8xWgf9qmEjjjfiX7P41r7x4yvH4ozJ6kWG2EhzxM5Embg9aKZtVLJdMAMCt1+/HXz5vEsfnc5geHeyeKwzTCRzxM5HGG/Ev2hH/ZEZDWldx5d7RPh4Zw2xdWPiZSCNHKxpmFYs5a7ThRJrr2xlmI7DwM5Em7o34c1bEP5HhVgYMsxFY+JlIo3si/oV8CZqqIKNzaophNgILPxNpZMRftiP+ybQGosFuqcswG4WFn4k09R7/RIb9fYbZKCz8TKTRvBF/3mB/n2G6AAs/E2lkHX/JtKweruhhmI3Dws9EGt3TpG0hV3LGGzIM0zks/EykkRH/SqGMklllq4dhugALPxNppMc/u7oOgBdvMUw32HThJ6KLiehhz781InrPZh8HszWQwn9utQiAF28xTDfY9JUwQoinAVwNAEQUA3AGwD9t9nEwW4O43ZZ51hb+SS7nZJgN02+r59UAjgshTvb5OJiIIuv4z67YVg9H/AyzYfot/LcC+EKfj4GJMESEeIwwbzdoG0+z8DPMRumb8BORBuC1AL4U8PN3EtFRIjo6Pz+/uQfHRAotpkAIYEhXIzVcnWG2Kv2M+H8WwINCiPN+PxRC3CaEOCKEODI1NbXJh8ZECZngZZuHYbpDP4X/zWCbhwmBrOXnPj0M0x36IvxElAJwI4A7+vH6zNbCifjZ32eYrtCXxuZCiAKAiX68NrP10DjiZ5iu0u+qHoZpiYz4uU8Pw3QHFn4m8rDVwzDdhYWfiTyc3GWY7sLCz0Qe1+PniJ9hugELPxN5XI+fI36G6QYs/Ezkcawe9vgZpiuw8DORR1cVKASMplj4GaYbsPAzkUdTFYylNMQU6vehMMy2oC8LuBimHW69bh9eeAGv92OYbsHCz0SeGy6YwA0s/AzTNdjqYRiGGTBY+BmGYQYMFn6GYZgBg4WfYRhmwGDhZxiGGTBY+BmGYQYMFn6GYZgBg4WfYRhmwCAhRL+PoSVENA/gZIiHTgJY6OAlRgCsbqPn8Hmw6OQ8dPI6nT6Pz8PmPmezzkOUzt0BIcRUw1YhxLb5B+Boh8+7bZs9h89Dh+ehk9fh88DnYSudOyEEWz02X9tmz+mUKL+nzToPnb4On4fOn7fdzkPUz93WsHrCQkRHhRBH+n0c/YbPgwWfBws+DxZ8Hly2W8R/W78PICLwebDg82DB58GCz4PNtor4GYZhmNZst4ifYRiGaUHkhZ+IPkVEc0T0mGfb1UR0HxE9TERHieh6e7tGRJ8momNE9AgRvcLznLuJ6Gn7OQ8T0Y7NfzedE3AeriKie+33+zUiGra330hED9jbHyCiV3me83x7+0+I6ONEtKXGWrVzHjw/309EOSL6Tc+2Qfo8vMXzPh8moioRXW3/bMt+Hto8B9tWGzqik1KgzfwH4GUArgXwmGfbvwH4WfvrmwDcbX/9bgCftr/eAeABAIr9/d0AjvT7/XT5PPwHgJfbX78dwB/aX18DYNr++jCAM57n3A/ghQAIwDfledwq/9o5D56ffxnAlwD8pmfbwHwe6p53BYBnt8Pnoc2/iW2rDZ38i3zEL4S4B8BS/WYAMqobATBrf30ZgLvs580BWAGwLbL4AefhYgD32F9/C8Ab7cc+JISQ5+RxAAki0oloN4BhIcS9wvrEfxbA63t+8F2knfMAAET0egDPwjoP24Z2z4OHNwP4AgBs9c9Dm+dg22pDJ0Re+AN4D4A/JaLTAD4K4Hft7Y8AeB0RqUR0CMDzAezzPO/T9q3cB7bSLW0THgPwWvvrN6H2vUreCOAhIUQJwB4AM56fzdjbtjq+54GI0gB+B8DvBzxvED8P/wW28GN7fh6CzsGgaUNTtqrw/zqA9woh9gF4L4C/tbd/CtaH9yiA/w/ADwGY9s/eIoS4AsBL7X+/tJkH3CPeDuDdRPQAgCEAhveHRHQ5gP8D4FflJp99bIeyrqDz8PsA/lwIkfN5ziB+Hm4AUBBCSE98O34egs7BoGlDU7bqsPW3Afgf9tdfAvBJABBCmLAuBAAAIvohgGfsn52x/88S0ecBXA/r1nbLIoR4CsDPAAARXQTgNfJnRLQXwD8BeKsQ4ri9eQbAXs8u9sK1ybYsTc7DDQBuIaI/ATAKoEpERSHEXw7a58HmVrjRPrANPw9B52DQtKEVWzXinwXwcvvrV8H+BRJRyr69BxHdCMAUQjxh395N2tvjAG6GdUu4pZHVB0SkAPg9AJ+wvx8F8A0AvyuE+IF8vBDiLIAsEb3Avp19K4A7N/u4u03QeRBCvFQIcVAIcRBWlPe/hRB/OWifB8+2NwH4R7ltO34emvxNDJQ2tCLyET8RfQHAKwBMEtEMgA8B+BUAf0FEKoAigHfaD98B4F+JqArgDNxbNt3eHgcQA/BtAH+zaW+iCwSchwwRvdt+yB0APm1//d8APA/AB4joA/a2n7GTWr8O4DMAkrCqOL65KW+gS7R5HoIYtM8DYFXAzAghnq3b1Zb9PLR5DratNnQCr9xlGIYZMLaq1cMwDMN0CAs/wzDMgMHCzzAMM2Cw8DMMwwwYLPwMwzADBgs/w9RBRKNE9C7762kiur3fx8Qw3YTLORmmDiI6CODrQojD/T4WhukFkV/AxTB94CMALiSih2GtCr9UCHGYiH4ZVvfKGKx21x8DoMFaDFQCcJMQYomILgTwfwFMASgA+BW7lQDDRAK2ehimkfcBOC6EuBrAb9X97DCAX4DVz+WPYTU9uwbAvbBaHgDWbNffEEI8H8BvAvirzThohgkLR/wM0x7fEUJkYfW4WQXwNXv7MQBXElEGwIsAfMnT3Vff/MNkmGBY+BmmPUqer6ue76uw/p4UACv23QLDRBK2ehimkSysXu5tI4RYA3CCiN4EAGRxVTcPjmE2Cgs/w9QhhFgE8AN7iPefdrCLtwB4BxE9Amvk4+u6eXwMs1G4nJNhGGbA4IifYRhmwGDhZxiGGTBY+BmGYQYMFn6GYZgBg4WfYRhmwGDhZxiGGTBY+BmGYQYMFn6GYZgB4/8HRPLjsf4TEdoAAAAASUVORK5CYII=\n", + "text/plain": [ + "<Figure size 432x288 with 1 Axes>" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "# Import Meteostat library and dependencies\n", + "from datetime import datetime\n", + "import matplotlib.pyplot as plt\n", + "from meteostat import Stations, Daily\n", + "import pandas as pd\n", + "\n", + "# Set coordinates of Hamburg Fuhlsbüttel\n", + "lat = 53.5510\n", + "lon = 9.9936\n", + "\n", + "# Set time period\n", + "start = datetime(1891, 1, 1)\n", + "end = datetime(2018, 12, 31)\n", + "\n", + "# Get closest weather station to Vancouver, BC\n", + "stations = Stations()\n", + "stations = stations.nearby(lat, lon)\n", + "stations = stations.inventory('daily', (start, end))\n", + "station = stations.fetch(1)\n", + "\n", + "# Get daily data for 2018 at the selected weather station\n", + "data = Daily(station, start, end)\n", + "data = data.fetch()\n", + "\n", + "# Calc Daily Mean\n", + "yearly_mean_df = data.groupby(pd.Grouper(freq='1Y')).mean()\n", + "\n", + "# Plot line chart including average, minimum and maximum temperature\n", + "yearly_mean_df.plot(y=['tmin', 'tmax'], ylabel = 'Degrees C', color=['green', 'red'])\n", + "plt.ylim(10., 15.)\n", + "plt.show()\n", + "\n", + "# Plot line chart including average, minimum and maximum temperature\n", + "yearly_mean_df.plot(y=['tavg'], ylabel = 'Degrees C')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "<div>\n", + "<style scoped>\n", + " .dataframe tbody tr th:only-of-type {\n", + " vertical-align: middle;\n", + " }\n", + "\n", + " .dataframe tbody tr th {\n", + " vertical-align: top;\n", + " }\n", + "\n", + " .dataframe thead th {\n", + " text-align: right;\n", + " }\n", + "</style>\n", + "<table border=\"1\" class=\"dataframe\">\n", + " <thead>\n", + " <tr style=\"text-align: right;\">\n", + " <th></th>\n", + " <th>name</th>\n", + " <th>country</th>\n", + " <th>region</th>\n", + " <th>wmo</th>\n", + " <th>icao</th>\n", + " <th>latitude</th>\n", + " <th>longitude</th>\n", + " <th>elevation</th>\n", + " <th>timezone</th>\n", + " <th>hourly_start</th>\n", + " <th>hourly_end</th>\n", + " <th>daily_start</th>\n", + " <th>daily_end</th>\n", + " <th>distance</th>\n", + " </tr>\n", + " <tr>\n", + " <th>id</th>\n", + " <th></th>\n", + " <th></th>\n", + " <th></th>\n", + " <th></th>\n", + " <th></th>\n", + " <th></th>\n", + " <th></th>\n", + " <th></th>\n", + " <th></th>\n", + " <th></th>\n", + " <th></th>\n", + " <th></th>\n", + " <th></th>\n", + " <th></th>\n", + " </tr>\n", + " </thead>\n", + " <tbody>\n", + " <tr>\n", + " <th>10147</th>\n", + " <td>Hamburg Airport</td>\n", + " <td>DE</td>\n", + " <td>HH</td>\n", + " <td>10147</td>\n", + " <td>EDDH</td>\n", + " <td>53.6333</td>\n", + " <td>10.0000</td>\n", + " <td>16.0</td>\n", + " <td>Europe/Berlin</td>\n", + " <td>1949-01-01</td>\n", + " <td>2020-12-15</td>\n", + " <td>1891-01-01</td>\n", + " <td>2020-12-03</td>\n", + " <td>9.161085e+03</td>\n", + " </tr>\n", + " <tr>\n", + " <th>10162</th>\n", + " <td>Schwerin</td>\n", + " <td>DE</td>\n", + " <td>MV</td>\n", + " <td>10162</td>\n", + " <td><NA></td>\n", + " <td>53.6500</td>\n", + " <td>11.3833</td>\n", + " <td>59.0</td>\n", + " <td>Europe/Berlin</td>\n", + " <td>1950-01-01</td>\n", + " <td>2020-12-15</td>\n", + " <td>1890-01-01</td>\n", + " <td>2020-12-04</td>\n", + " <td>9.235691e+04</td>\n", + " </tr>\n", + " <tr>\n", + " <th>10224</th>\n", + " <td>Bremen</td>\n", + " <td>DE</td>\n", + " <td>HB</td>\n", + " <td>10224</td>\n", + " <td>EDDW</td>\n", + " <td>53.0500</td>\n", + " <td>8.8000</td>\n", + " <td>3.0</td>\n", + " <td>Europe/Berlin</td>\n", + " <td>1926-01-01</td>\n", + " <td>2020-12-15</td>\n", + " <td>1890-01-01</td>\n", + " <td>2020-12-03</td>\n", + " <td>9.692615e+04</td>\n", + " </tr>\n", + " <tr>\n", + " <th>D2578</th>\n", + " <td>Kirchdorf/Poel</td>\n", + " <td>DE</td>\n", + " <td>MV</td>\n", + " <td><NA></td>\n", + " <td><NA></td>\n", + " <td>53.9995</td>\n", + " <td>11.4341</td>\n", + " <td>12.0</td>\n", + " <td>Europe/Berlin</td>\n", + " <td>2004-07-01</td>\n", + " <td>2020-12-02</td>\n", + " <td>1865-01-01</td>\n", + " <td>2020-12-04</td>\n", + " <td>1.069908e+05</td>\n", + " </tr>\n", + " <tr>\n", + " <th>10361</th>\n", + " <td>Magdeburg</td>\n", + " <td>DE</td>\n", + " <td>ST</td>\n", + " <td>10361</td>\n", + " <td>EDBM</td>\n", + " <td>52.1333</td>\n", + " <td>11.6000</td>\n", + " <td>79.0</td>\n", + " <td>Europe/Berlin</td>\n", + " <td>1951-01-01</td>\n", + " <td>2020-12-15</td>\n", + " <td>1881-01-01</td>\n", + " <td>2020-12-04</td>\n", + " <td>1.910266e+05</td>\n", + " </tr>\n", + " <tr>\n", + " <th>...</th>\n", + " <td>...</td>\n", + " <td>...</td>\n", + " <td>...</td>\n", + " <td>...</td>\n", + " <td>...</td>\n", + " <td>...</td>\n", + " <td>...</td>\n", + " <td>...</td>\n", + " <td>...</td>\n", + " <td>...</td>\n", + " <td>...</td>\n", + " <td>...</td>\n", + " <td>...</td>\n", + " <td>...</td>\n", + " </tr>\n", + " <tr>\n", + " <th>94964</th>\n", + " <td>Bushy Park</td>\n", + " <td>AU</td>\n", + " <td>TAS</td>\n", + " <td>94964</td>\n", + " <td><NA></td>\n", + " <td>-42.7000</td>\n", + " <td>146.8833</td>\n", + " <td>35.0</td>\n", + " <td>Australia/Hobart</td>\n", + " <td>NaT</td>\n", + " <td>NaT</td>\n", + " <td>1889-04-03</td>\n", + " <td>2020-12-04</td>\n", + " <td>1.855174e+07</td>\n", + " </tr>\n", + " <tr>\n", + " <th>95974</th>\n", + " <td>Fingal</td>\n", + " <td>AU</td>\n", + " <td>TAS</td>\n", + " <td>95974</td>\n", + " <td><NA></td>\n", + " <td>-41.6333</td>\n", + " <td>147.9667</td>\n", + " <td>233.0</td>\n", + " <td>Australia/Hobart</td>\n", + " <td>NaT</td>\n", + " <td>NaT</td>\n", + " <td>1888-10-03</td>\n", + " <td>2020-11-29</td>\n", + " <td>1.857037e+07</td>\n", + " </tr>\n", + " <tr>\n", + " <th>94979</th>\n", + " <td>Lake Leake</td>\n", + " <td>AU</td>\n", + " <td>TAS</td>\n", + " <td>94979</td>\n", + " <td><NA></td>\n", + " <td>-42.0000</td>\n", + " <td>147.7833</td>\n", + " <td>580.0</td>\n", + " <td>Australia/Hobart</td>\n", + " <td>NaT</td>\n", + " <td>NaT</td>\n", + " <td>1890-01-05</td>\n", + " <td>2020-12-04</td>\n", + " <td>1.858111e+07</td>\n", + " </tr>\n", + " <tr>\n", + " <th>94970</th>\n", + " <td>Hobart Regional Office</td>\n", + " <td>AU</td>\n", + " <td>TAS</td>\n", + " <td>94970</td>\n", + " <td><NA></td>\n", + " <td>-42.8833</td>\n", + " <td>147.3167</td>\n", + " <td>51.0</td>\n", + " <td>Australia/Hobart</td>\n", + " <td>NaT</td>\n", + " <td>NaT</td>\n", + " <td>1882-04-01</td>\n", + " <td>2020-12-04</td>\n", + " <td>1.860456e+07</td>\n", + " </tr>\n", + " <tr>\n", + " <th>94967</th>\n", + " <td>Cape Bruny</td>\n", + " <td>AU</td>\n", + " <td>TAS</td>\n", + " <td>94967</td>\n", + " <td><NA></td>\n", + " <td>-43.5000</td>\n", + " <td>147.1500</td>\n", + " <td>55.0</td>\n", + " <td>Australia/Hobart</td>\n", + " <td>NaT</td>\n", + " <td>NaT</td>\n", + " <td>1880-01-03</td>\n", + " <td>2020-12-04</td>\n", + " <td>1.863516e+07</td>\n", + " </tr>\n", + " </tbody>\n", + "</table>\n", + "<p>229 rows × 14 columns</p>\n", + "</div>" + ], + "text/plain": [ + " name country region wmo icao latitude \\\n", + "id \n", + "10147 Hamburg Airport DE HH 10147 EDDH 53.6333 \n", + "10162 Schwerin DE MV 10162 <NA> 53.6500 \n", + "10224 Bremen DE HB 10224 EDDW 53.0500 \n", + "D2578 Kirchdorf/Poel DE MV <NA> <NA> 53.9995 \n", + "10361 Magdeburg DE ST 10361 EDBM 52.1333 \n", + "... ... ... ... ... ... ... \n", + "94964 Bushy Park AU TAS 94964 <NA> -42.7000 \n", + "95974 Fingal AU TAS 95974 <NA> -41.6333 \n", + "94979 Lake Leake AU TAS 94979 <NA> -42.0000 \n", + "94970 Hobart Regional Office AU TAS 94970 <NA> -42.8833 \n", + "94967 Cape Bruny AU TAS 94967 <NA> -43.5000 \n", + "\n", + " longitude elevation timezone hourly_start hourly_end \\\n", + "id \n", + "10147 10.0000 16.0 Europe/Berlin 1949-01-01 2020-12-15 \n", + "10162 11.3833 59.0 Europe/Berlin 1950-01-01 2020-12-15 \n", + "10224 8.8000 3.0 Europe/Berlin 1926-01-01 2020-12-15 \n", + "D2578 11.4341 12.0 Europe/Berlin 2004-07-01 2020-12-02 \n", + "10361 11.6000 79.0 Europe/Berlin 1951-01-01 2020-12-15 \n", + "... ... ... ... ... ... \n", + "94964 146.8833 35.0 Australia/Hobart NaT NaT \n", + "95974 147.9667 233.0 Australia/Hobart NaT NaT \n", + "94979 147.7833 580.0 Australia/Hobart NaT NaT \n", + "94970 147.3167 51.0 Australia/Hobart NaT NaT \n", + "94967 147.1500 55.0 Australia/Hobart NaT NaT \n", + "\n", + " daily_start daily_end distance \n", + "id \n", + "10147 1891-01-01 2020-12-03 9.161085e+03 \n", + "10162 1890-01-01 2020-12-04 9.235691e+04 \n", + "10224 1890-01-01 2020-12-03 9.692615e+04 \n", + "D2578 1865-01-01 2020-12-04 1.069908e+05 \n", + "10361 1881-01-01 2020-12-04 1.910266e+05 \n", + "... ... ... ... \n", + "94964 1889-04-03 2020-12-04 1.855174e+07 \n", + "95974 1888-10-03 2020-11-29 1.857037e+07 \n", + "94979 1890-01-05 2020-12-04 1.858111e+07 \n", + "94970 1882-04-01 2020-12-04 1.860456e+07 \n", + "94967 1880-01-03 2020-12-04 1.863516e+07 \n", + "\n", + "[229 rows x 14 columns]" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "stations.fetch(0)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "test_env", + "language": "python", + "name": "test_env" + }, + "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.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/notebooks/use-case_count_summer_days_cmip6.ipynb b/notebooks/use-case_count_summer_days_cmip6.ipynb index 0b6b25a148967d303f32ede0a82e8b2aa1b20ef1..60a10c8a6a5987a6f6101e45f109b6f616461bc9 100644 --- a/notebooks/use-case_count_summer_days_cmip6.ipynb +++ b/notebooks/use-case_count_summer_days_cmip6.ipynb @@ -197,7 +197,7 @@ "source": [ "# Store the name of the model we chose in a variable named \"climate_model\"\n", "\n", - "climate_model = \"MPI-ESM1-2-LR\" # here we choose Max-Plack Institute's Earth Sytem Model in high resolution\n", + "climate_model = \"MPI-ESM1-2-HR\" # here we choose Max-Plack Institute's Earth Sytem Model in high resolution\n", "\n", "# This is how we tell intake what data we want\n", "\n", @@ -258,7 +258,7 @@ "outputs": [], "source": [ "# Select the file that contains the year we selected in the drop down menu above, e.g. 2015\n", - "selected_file = query_result_df_m[(year_box.value >= query_result_df[\"start_year\"]) & (\n", + "selected_file = query_result_df[(year_box.value >= query_result_df[\"start_year\"]) & (\n", " year_box.value <= query_result_df[\"end_year\"])]\n", "\n", "# Path of the file that contains the selected year \n", diff --git a/notebooks/use-case_dask_for_climate_data.ipynb b/notebooks/use-case_dask_for_climate_data.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..e66fa9a0fdda6cee3b5fff0d15dc73f7ca16e23d --- /dev/null +++ b/notebooks/use-case_dask_for_climate_data.ipynb @@ -0,0 +1,359 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# How to use Dask for Climate Data Processing?\n", + "\n", + "This tutorial builds requires the skills, which have learned in the summer days tutorial (provide link)\n", + "\n", + "1. What is Dask?\n", + " 1. Parallelism\n", + " 1. Overview\n", + " 1. `dask.delayed`\n", + "2. Process climate data with dask?\n", + "3. Common mistakes" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. What is Dask?\n", + "Dask is an open source library for parallel computing written in Python. It is used to process larger-than memory datasets (e.g. large climate data sets). All information can be found here: https://docs.dask.org" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1C. Parallelism\n", + "- use Maria's metaphor" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1B. Dask Overview" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### 1C. `dask.delayed`" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let us start with an easy example" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from dask.distributed import Client\n", + "\n", + "client = Client(n_workers=4)\n", + "\n", + "client" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from dask import delayed\n", + "import time" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Not Parallel" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "%%time\n", + "\n", + "def inc(x):\n", + " time.sleep(0.5)\n", + " return x + 1\n", + "\n", + "def double(x):\n", + " time.sleep(0.5)\n", + " return 2 * x\n", + "\n", + "def add(x, y):\n", + " time.sleep(0.5)\n", + " return x + y\n", + "\n", + "data = list(range(4))\n", + "\n", + "output = []\n", + "for x in data:\n", + " a = inc(x)\n", + " b = double(x)\n", + " c = add(a, b)\n", + " output.append(c)\n", + "\n", + "total = sum(output)\n", + "total" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Parallel" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "@delayed\n", + "def inc(x):\n", + " time.sleep(0.5)\n", + " return x + 1\n", + "\n", + "@delayed\n", + "def double(x):\n", + " time.sleep(0.5)\n", + " return 2 * x\n", + "\n", + "#@delayed\n", + "def add(x, y):\n", + " time.sleep(0.5)\n", + " return x + y\n", + "\n", + "data = list(range(4))\n", + "\n", + "output = []\n", + "for x in data:\n", + " a = inc(x)\n", + " b = double(x)\n", + " c = add(a, b)\n", + " output.append(c)\n", + "\n", + "total_delayed = sum(output) #also delay sum because it is a function\n", + "#%time total_delayed.compute()\n", + "total_delayed.visualize()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "total_delayed.visualize()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "client.close()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Looks good in theory, right? Now, let us apply our knowldege to climate model data." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 2. Process Climate Data with Dask\n", + "- load with intake\n", + "- think about processing\n", + "- compare conventional with dask" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import intake\n", + "import xarray as xr\n", + "#from import Basemap, cm\n", + "import matplotlib.pyplot as plt\n", + "from netCDF4 import Dataset as open_ncfile\n", + "import numpy as np" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Path to catalog descriptor on the DKRZ server\n", + "col_url = \"/work/ik1017/Catalogs/mistral-cmip6.json\"\n", + "\n", + "# Open the catalog with the intake package and name it \"col\" as short for \"collection\"\n", + "col = intake.open_esm_datastore(col_url)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Store the name of the model we chose in a variable named \"climate_model\"\n", + "\n", + "climate_model = \"MPI-ESM1-2-HR\" # here we choose Max-Plack Institute's Earth Sytem Model in high resolution\n", + "\n", + "# This is how we tell intake what data we want\n", + "\n", + "query = dict(\n", + " source_id = \"MPI-ESM1-2-HR\", # the model \n", + " variable_id = \"tasmax\", # temperature at surface, maximum\n", + " table_id = \"day\", # daily maximum\n", + " experiment_id = \"historical\", # what we selected in the drop down menu,e.g. SSP2.4-5 2015-2100\n", + " member_id = \"r10i1p1f1\", # \"r\" realization, \"i\" initialization, \"p\" physics, \"f\" forcing\n", + " time_range =\"20100101-20141231\",\n", + ")\n", + "\n", + "# Intake looks for the query we just defined in the catalog of the CMIP6 data pool at DKRZ\n", + "cat = col.search(**query)\n", + "\n", + "# Show query results\n", + "cat.df" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ds = xr.open_dataset(cat.df['path'][0])\n", + "ds" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "var = ds.variables['tasmax'][0,:,:]\n", + "lat = ds.variables['lat'][:]\n", + "lon = ds.variables['lon'][:]\n", + "\n", + "#-- create figure and axes instadses\n", + "dpi = 100\n", + "fig = plt.figure(figsize=(1100/dpi, 1100/dpi), dpi=dpi)\n", + "ax = fig.add_axes([0.1,0.1,0.8,0.9])\n", + "\n", + "#-- create map\n", + "map = Basemap(projection='cyl',llcrnrlat= -90.,urcrnrlat= 90.,\\\n", + " resolution='c', llcrnrlon=-180.,urcrnrlon=180.)\n", + "\n", + "#-- draw coastlines, state and country boundaries, edge of map\n", + "map.drawcoastlines()\n", + "map.drawstates()\n", + "map.drawcountries()\n", + "\n", + "#-- create and draw meridians and parallels grid lines\n", + "map.drawparallels(np.arange( -90., 90.,30.),labels=[1,0,0,0],fontsize=10)\n", + "map.drawmeridians(np.arange(-180.,180.,30.),labels=[0,0,0,1],fontsize=10)\n", + "\n", + "#-- convert latitude/longitude values to plot x/y values\n", + "x, y = map(*np.meshgrid(lon,lat))\n", + "\n", + "#-- contour levels\n", + "clevs = np.arange(210,320,5)\n", + "\n", + "#-- draw filled contours\n", + "cnplot = map.contourf(x,y,var,clevs,cmap=plt.cm.jet)\n", + "\n", + "#-- add colorbar\n", + "cbar = map.colorbar(cnplot,location='bottom',pad=\"10%\") #-- pad: distadse between map and colorbar\n", + "cbar.set_label('deg K') #-- add colorbar title string\n", + "\n", + "#-- add plot title\n", + "plt.title('Temperature')\n", + "\n", + "#-- displa on screen\n", + "#plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# <font color='darkred'>Now use Gradient Operator in parallel</font>" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# <font color='darkgreen'>Take home message</font>\n", + "\n", + "Parallelism brings extra complexiity and often it is not necessary for your problems. Before using Dask you may want try alternatives:\n", + "- use better algorithms or data structures\n", + "- better file formats\n", + "- compiled code\n", + "- sampling\n", + "- profile your code" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 unstable (using the module python3/unstable)", + "language": "python", + "name": "python3_unstable" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.8" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/notebooks/use-case_plot-unstructured_psyplot_cmip6.ipynb b/notebooks/use-case_plot-unstructured_psyplot_cmip6.ipynb new file mode 100644 index 0000000000000000000000000000000000000000..3904b570cf8b2bed96f0e1994834a83bc462bd2c --- /dev/null +++ b/notebooks/use-case_plot-unstructured_psyplot_cmip6.ipynb @@ -0,0 +1,229 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Plot Earth System Model data on unstructured grids with psyplot\n", + "\n", + "This notebook introduces you to the `mapplot` function of the package `psyplot`.\n", + "It is suitable to plot maps from data on unstructured grids like the ones from ICON and FESOM.\n", + "\n", + "We therefore search for the corresponding data in the CMIP6 data pool with intake-esm.\n", + "Afterwards, we open a file with `xarray` and configure the opened xarray dataset as well as psyplot for a map plot.\n", + "\n", + "In the end, we discuss the functions of `psyplot.project.plot.mapplot`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import psyplot.project as psy\n", + "import matplotlib as mpl\n", + "import xarray as xr\n", + "import intake" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We open a swift catalog from dkrz cloud which is accessible without authentication." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "col_url = \"https://swift.dkrz.de/v1/dkrz_a44962e3ba914c309a7421573a6949a6/intake-esm/mistral-cmip6.json\"\n", + "col = intake.open_esm_datastore(col_url)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this example, we aim at plotting the Sea Surface Temperature of the upper boundary of the liquid ocean, including temperatures below sea-ice and floating ice shelves from AWI.\n", + "We therefore search for `tos` in the catalog for monthly frequency. We use 1 realization of 1 experiment only." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "tos=col.search(source_id=\"AWI-CM-1-1-MR\",\n", + " experiment_id=\"ssp370\",\n", + " variable_id=\"tos\",\n", + " table_id=\"Omon\",\n", + " member_id=\"r1i1p1f1\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We now open the file on the mistral file system." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dset = xr.open_dataset(tos.df[\"path\"].to_list()[0])\n", + "#dset = xr.open_mfdataset(tos.df[\"path\"].to_list())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In order to make `tos` plottable, we set the following configuration.\n", + "- The `CDI_grid_type` is a keyword for `psyplot`.\n", + "- Coordinates are not fully recognized by `xarray` so that we have to add some manually (version from Dec 2020)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dset[\"tos\"][\"CDI_grid_type\"]=\"unstructured\"\n", + "coordlist=[\"vertices_latitude\", \"vertices_longitude\", \"lat_bnds\", \"lon_bnds\"]\n", + "dset=dset.set_coords([coord for coord in dset.data_vars if coord in coordlist])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This is based on the example from:\n", + "https://psyplot.readthedocs.io/projects/psy-maps/en/latest/examples/example_ugrid.html#gallery-examples-example-ugrid-ipynb" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "psy.rcParams['plotter.maps.xgrid'] = False\n", + "psy.rcParams['plotter.maps.ygrid'] = False\n", + "mpl.rcParams['figure.figsize'] = [10, 8.]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "iconplot11=psy.plot.mapplot(\n", + " dset, name=\"tos\", cmap='rainbow',\n", + " clabel=dset[\"tos\"].description,\n", + " stock_img=True, lsm='50m')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We now do the same with a smaller subset to see the fine resolution of the AWI ocean model FESOM.\n", + "The subsetting is required because the plotting takes too long otherwise." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dset2 = dset.isel(time=slice(1,3)).where( (dset.lon > -10. ) &\n", + " (dset.lon < 50. ) &\n", + " (dset.lat > 40. ) &\n", + " (dset.lat < 70. ), drop=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dset2.to_netcdf(\"/home/dkrz/k204210/test.nc\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dset2=xr.open_dataset(\"/home/dkrz/k204210/test.nc\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dset2[\"tos\"][\"CDI_grid_type\"]=\"unstructured\"\n", + "coordlist=[\"vertices_latitude\", \"vertices_longitude\", \"lat_bnds\", \"lon_bnds\"]\n", + "dset2=dset2.set_coords([coord for coord in dset2.data_vars if coord in coordlist])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "iconplot12=psy.plot.mapplot(\n", + " dset2, name=\"tos\", cmap='rainbow',\n", + " lonlatbox='Ireland',\n", + " clabel=dset[\"tos\"].description,\n", + " stock_img=True,\n", + " lsm='50m',\n", + " datagrid=dict(c='b', lw=0.2))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 unstable (using the module python3/unstable)", + "language": "python", + "name": "python3_unstable" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.8" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +}