Skip to content

Instantly share code, notes, and snippets.

@kmuehlbauer
Last active January 13, 2024 14:18
Show Gist options
  • Save kmuehlbauer/645e42a53b30752230c08c20a9c964f9 to your computer and use it in GitHub Desktop.
Save kmuehlbauer/645e42a53b30752230c08c20a9c964f9 to your computer and use it in GitHub Desktop.
xarray and pysteps, just good friends
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# HowTo use xarray withon pysteps\n",
"\n",
"This just shows how xarray data model could be used to simplify reading/processing/plotting."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Pysteps configuration file found at: /home/kai/miniconda3/envs/wradlib/lib/python3.7/site-packages/pysteps/pystepsrc\n",
"\n"
]
}
],
"source": [
"import xarray as xr\n",
"import cartopy.crs as ccrs\n",
"import matplotlib.pyplot as pl\n",
"import pysteps"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"filename = '/home/data/kai/projects/pysteps-data/radar/bom/prcp-cscn/2/2018/06/16/2_20180616_100000.prcp-cscn.nc'"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## pysteps original"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
" R, _, metadata = pysteps.io.import_bom_rf3(filename)"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<cartopy.mpl.geoaxes.GeoAxesSubplot at 0x7f71d035b780>"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAASwAAADnCAYAAACkJWu2AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnXecVNX1wL9nG0tvS+9NEKS3FRsCVoINjSUaW2K6JVGjMWKCyS+amChRE5PYE6NCsCEKtig2QPoCSkAEpEgTlu1l5vz+eDO7b2bezLypO495Xz7Lztx37313Z+4799x7zz1HVBUXFxcXJ5DT1A1wcXFxsYsrsFxcXByDK7BcXFwcgyuwXFxcHIMrsFxcXBxDXqSLvXv31C+/3JWutri4ZCvbVbVvUzfCCUgkswYR0UOHNgekHThwhPLyavr27ZzqtiWVzZt306VLO9q0adHUTYmJlSs/Z+zYAU3djJioqqpl69avGDasd1M3JSZ27/4agO7dO6T1vu3bD0JVJa03dShxTAmd+rkKrslZOnFmPxFnNjtriFlgiYATjU3djugSDaNbux0lk4lTYKWiKanHiYLWmahDBwintjt7iHNK6D74LuFx6rjgaliZTxZpWG5HdImOq2FlNnEILGcuXjtX0LqkC3fJIPOJ03DUqV+sU9vtki7EVbEymoiGo1a4GpaLHZz44GeDhiUilwL9gXygXFXva+ImxURca1jO1FSc9wA5FSc/9w6Us7FyJrCoqRsRL3HtEjq5Q7qkA2d2kEzcJRSRviLymYg8KiLrReQZEZkmIh+KyGYRmSAivxKRp0TkDRHZJiIXiMjvRaRERBaJSL6vLgFGAat81Q8VkXdFZKuIXN9kf2QMZJXhqBPb7ZJOMtYOayAwBxgBDAEuA04EbgZ+4cszAJgOnAv8C/ivqg4HqnzpAKOBtdr4IAwBzgAmAHf5BVsmE9cuoYtLNJzYTTJ4PPtCVUtU1QtsAN72CZ0SoK8vz+uqWudLy6Vx2mfOcybwuqnehapao6oHgH1Al5T+FUkgizQsBz5BDsWB3aOBDO0nNabXXtN7L40bZzUAPqFWZ9KizHlOB94IU6+HODbh0k1cZg1O7ZBObbczycgHPyJOHIjtIiJtgTxVPdjUbUmELJsSHr0d0iU5OLt/R+Q04K2mbkSixGGH5cyR6CjuiBlIxi5eRyQTu7WqbgOOM72/Ktw1U3or0+tfAYjIo8Cjwemm9yH1ZCJxzFmda9bg1Ha7pAtnClo7qOp3mroNySCqwPJ6vag2alUejwev10NtbZ0vHUAD8jT+Dnztz9eY35xmXd76upFunSfwvb9caWkllZU1VFXVBpUN1BhD6264EiLwzHlD0wOvRRKWkTTW6uo6Nm0KdFOdiOCN9YEMp5kGJktAem1tPRUV1WzZsgeQkHsa78WyHvN763yNeayum9vrf+3PF5wnOL2yspbc3Crf9yERyzW+tk4zlw18HZjH1fxjI6rAWrNmW8CH6/UqVVU1bNq0G7D6ckI7g3UHCewI0TpTcLp1Hn9d/k7QeC0vL4f8/FwKC/MtO7/Ve/O9Gt+HKxuaP7RMcNtDCU4vLa2kR4+O1pmTgpVwtshlIYStyymVlTV4PF6KitoE5bMS4MGDi1W+0DYGDl6B9ZgHvvCDX2iZ2tp6KitrGgZp64E13H0Ieh08iFvlDx0EXSITVWCNGdM/4H1NTR3/+99uhg/vk7JGpYKamjoKCwvo3LltUzclJnJzc2jVqrCpmxETOTk5HDlSRbt2LZu6KTFRX++hqKiN49odjIi0VdXSpm5HKsiaXUKnttupOPXzdmq7/YjIwBw4LCKOWESPlazZJXRxicbR0K3HwuYOwF7Dwt3Z0teCLDMcdWjDHYZzP2Zn7xKKyMBDGAZXVcb7o07LissflhNxaLMdijMllqpz+zcY2tUYDLXqNKD5UahlxalhObFDOtd+zIk48bl3Zr828GtXg3zv+3F0allZ49PdwLENd0kTThS0YGhXUwlUp04DRhhaVlRE5HER2Sci64PSfyIim0Rkg4j83pR+u4hs8V07Iyl/hA2yyFtDU7cge8hER3h2cGq7g7UrPzFqWU9iuJ8x13sqhn+tEao6DLjPlz4UuAQY5ivzFxHJTeRvsEvWmDWAkxeDXdKDMxfdrbQrP3a1LFVdAnwdlPwD4B5V9buu2edLPxd4zudL6wtgC4YTwJSTNWtYTha0TsSJH7cDuzUikreLUO3KTz8MKeSb2q0w/Vxno/pjgJNEZJmIvCci433pPYAvTfl2+tJSTtbsEoIzBa1LenFg/5Zu+XBxUfgMLx6EnbX8U1UfjLHuPKA9UAyMB+aKSH+slbm0PFxxxiV0cYmEMwcGd0ALYSfwghosx/BeWuRL72XK1xPYnY4GZY3Act7A6Vyc/Ny7/SSAl4ApACJyDFAAHABeAS4RkWYi0g9jRro8HQ3KeB/OyUPcEdQlIk7dJUwGIvIsMBkoEpGdwF3A48DjPlOHWuBKn6/4DSIyF9gI1AM/UlVPOtqZRQILnDpVcSZOfPCzt3+o6qVhLl0eJv9vgd+mrkXWuFNCl5TgxM/bOJrT1K1wiUTWCCz3aE76cPLn7MBdwqwiawSW2w/TiTMlVjavcVodzRGRDiLypohs9v1ub7rmjKM5TiabO6SLPbJYw3qSoKM5wG0YUaYHAW/73jvraI6Ly9FKNo9nYY7mnAs85Xv9FHCeKb1JjuZkzS5hFo+cace5i9fOPEtY0Bl6XxP+euE/gW1cKSJXmJL/rqp/j1J1F1XdA6Cqe0Sksy+9B7DUlC9zj+Y4mWweQdOP8578o9wO6ylVnZOkutyjOenBlVgu4VF1poaVQvaKSDcA32+/twb3aE6qMfx4NXUrXFwcxSvAlb7XVwIvm9LdozmpxR0604czRwan+3RPhDBHc+7B8NBwLbADuAhAVd2jOakmS/thk+AuujuPCEdzpobJ3yRHc5pMYIUL8238BuuQ4Ea6dZ7A98HXKitrqK6upbS0MuS6+R6NdYXabQWHUQ/OHynsezQbMOvLSl1dPXv2HAp4kBKZ2kZ+IEMvhstvpYn4k8rKqqmqqmH//iMW5SUgzXgtQe+t8oXmCUwTm9ckJI8/3etVPB5vgMC1qsel6YgqsFat2hrysFVW1rBy5ecJ39zcYRpfR+5cka5ZdXb/terqWurqPOzbVxr0EIV/gCLlC84bmt+6TGj5wDqCUQWPx5vykT+yENSw162EeE1NHXV1HiorayzLhg481vcJrtuqXOhAEzrQRbpmTq+qqmXDhi8brlsNfskkkwSgiLQDHgWOw/iQrwHOxrC58mIsuF+lqmlZXA9HVIE1alS/AIECsHLl54wZ0z+jPvBoHDhQRnl5FX37do6eOYPYu/cwPXp0cNRnffBgGUeOVNGnT6embkpMrFr1OSNH9k3LZ22eYWQIc4BFqnqhiBQALYANqnongIhcD8wCvt+EbYwusHJy7E8TMhmn7hI68bN26hpWOvtH8IyiKRGRNsDJwFUAqlqL4f/KTEsyYDclTrMGp3o+cGSjHYhzP+dMECAp4soIQSj6A/uBJ0RktYg8KiItAUTktyLyJfAtDA2rSYlLYBnfqXM7pYvLUUUzYECEn+aAYek+zvRjPpaTB4wB/qqqo4EKfAedVfUOVe0FPAP8OD1/UHiyRsNybsTqo3bEd8kcdgI7VXWZ7/1/MASYmX8DM9PaKgvi1rAyaLHQxcUlAVT1K+BLERnsS5oKbBQRc7jDc4DP0t64IOKyw3LiPN+dxqYXB3aRbOcnwDO+HcKtwNXAoz4h5gW208Q7hJCA4agTNSwHNtmRHOVeD45KVHUNMC4oOelTQBHpoKrBfrdsE+eU0InrQe4D5OKSASwTkXkicrbEMVXLml1Cp9phORP3g3YiIpLrM2t41fc+rE/3BDgG+DtwBbBFRP7PF6TVFnG7l3Hmw+/IRjsSdw3LkdwAfGp6b+nTPRF8Ye/f9B22/g6G25rlIvKeiBwfrXzWTAndByh9OK1vuICI9ASmY5wn9BPOp3si9+koIjeIyArgZozF/iLgZximExGJc5cQnKetOE/I+nHqUReXjCKaT/cHgFuB1qa0cD7dE+Fj4J/Aeaq605S+QkQeiVY4zl1C5z78LunClbBpozkwJML1lkAEn+4i8g1gn6quFJHJSW9fIL9U1blB979IVeep6r3RCmeN4agztULn4mqEjuIE4BwR2QY8B0wRkX8R3qd7Ilitg91ut3DWGI6Ck9dWFCdpLE4bzLIdVb0dn9DwaVg3q+rlIvIHjEXxewj06R4zInIWhn+tHiLyZ9OlNhhulm2RNYajThWyLi5NiKVP9zjZDazAOOKz0pReBtxkt5K4NSyHySsXFxcbqOq7wLu+1wcJ49M9jnrXAmtF5BlVta1RBZNFu4TO0wpdXI4WRGSuqn4TWC0i5gdRMMyzRtipJ4EpYbwlmwanzgjddrscJdzg+/2NRCpJwHDUYRLLNcVIG+7hZ2chIr1E5L8i8qmIbBCRG3zpvxKRXSKyxvdzdrz38NtzAQeAL1V1O4brwZHEEDU67imhMx9+RzbaxSXV1AM/U9VVItIaWCkib/qu3a+q9yXxXkuAk3znEt/GWIi/GMMFc1Ti1rCchgOb7GDcgcFJqOoeVV3le12GcZ6wR4puJ6paCVwAPKiq5wND7RbOGrMGcKpW6EzcASKNtAImRrjeBoh+NAcAEekLjAaWYRiU/lhEvo2hCf1MVQ8l2FrxHXL+FnCtL822HMqiw8/uE5QunNY3soRIQSgAEJFWwHzgRlU9AvwVI4zFKGAP8McktOMGDCPVF1V1g4j0B/5rt3ACoeqd1yudqBW6mwUu6UBE8jGE1TOq+gKAqu41Xf8H8Gqi91HVJRjrWP73W4Hr7ZbPKsNRI/x4YMTdwNDokcOi+65EDL0emBZYzjo9chmPx0tlZY1lQNvk4w/uGSWXNOYNTDOoq/Pg9So1NXVhywXeQwLeW+cLzpPcz8OZg1ly8Hn+fAz4VFX/ZErvZtrdOx9Yn4R7HYPhVqYvJvmjqlPslI8qsDZv3oP/IfV/qWVlVRw5UsnevYd9QoCQPNHSrB78VOL1KrW1daxe/UVDWrgHyPwwmB+Y4HzBeUOvW9UTSqSHr7a2ji+/PICIpGRdyM53YPUwBwtxc3p1dR2glJdXRShrZ8Aw8oUbMA7vOwJAu85tov8RNqisrGHlys9D0oP7SXDUZvN30/hagvJa5U/NdxoHJ2B4AC0RkTW+tF8Al4rIKIwvaxvwvSTcax7wCIbfLU+shaMKrM6d24Z8Abt2fU2rVoW0b98y6IP357GXBulbW6qqqmXr1r0MG9YrLfdLFqtXf8HgwT3IzY3bOWza2bPna1She/cOKbvHwfZGv+l4KDmjnterrFnzBWPG9A9IDxSijcLSSlM3pwUO2ObX1gN8U6KqH2BtOPdaCm5Xr6p/jbdwVIHVtm2LkLRmzfJo1iyPFi2axXvfJqLpO0csHGwvfPF/XzJiRJ+mbkpMpOMZLKEEgMkcl6Qa1VLbCdS2M0MdcjgLROSHwItAjT/RbiSdrFnDcpqx6/z7F7OUNVxW/M2mbkpGMvlQsgSV8VkjQv9TBkXPfJQiIo9jHJvZp6rHmdJ/ghGivh5YqKq3JnirK32/bzGlKdDfIm8ICcwzHPT0O5CZN53BHw79PFPWOI465t+/ODChYdkia3kSONOcICKnYvh1H6Gqw4CELd5VtZ/Fjy1hBQl5HI2nZFOS3b0x3Tjh4b+lveGRd+ZNZyDA/p1xx/d0PD5zg+AP4AfAPapa48uTsMdREWkhIr8Ukb/73g/yuWi2hXs0x0G0f/bN6JkygEw//Dz//sUsnb2G4lmjGjStc3409WifNFwpIitMP9fZKHMMxrm/Zb4wXOOT0I4ngFpgku/9TuA3dgtn2dEc57XZz/z7FzP61XFwaVO3xPnMvOkMZt50BtA4NVRVOvfuyPz7FzdccwrlGOdowmEYf4QPQhGBPKA9UAyMx/A+2l8Te5AGqOrFInIpgKpWxRIBOouO5jR1CxJBWTp7De+83bGpG2ITZ3UOv9ASwXHCKsXsBF7wBT9dDngxYggmQq2INMfXSURkAKbdwmhkTah6px9x+cOhnzd1E2IiUweIkMV2DCFlxH7M0EY3HS8BU6DBQr0Aw59VIvwKWAT0EpFnMFzM2O7cce4SOvvhdzGYf/9iDrYX2t8ysKmbkjbMGlSg8MruDi0iz2IEOR0sIjt9gSceB/qLyHqM8F9XJjgdRFXfwHAtcxXwLDBOVVN7+NmNS5heUvVRz7zpDLhJOcSWpNab6YvufszCK9s1LFUNtzp6eTLvIyJvq+pUYKFFWlSyyuOoE9vs4nI0ICKFQAugyOdt1D86tAG6260ngSmh057+7B09XaKjan00J1sQkcdFZJ9v+udP6yAib4rIZt/v9gnc4nsY8QiH+H77f14GHrZbibtL6JISnPZ5O2Uam0KeJMjSHSOs/NuqOghjcdwqzLwtVHWOqvbDiCrd32TlPlJVH7JbjxuX0CUFOPNzdpqQTSaqusTnHtnMucBk3+unMAKsJrRdraoPisgkQv1hPW2nvBuX0CXpOK1vgDuYhaGL34Gfqu4Rkc6JVigi/8Rwu7yGRn9YCqROYLlxCV2ORpy6S1gBrIpwvcz4ZSsIRRoYBwyN1zwiq3YJnYhTHyKntTsL+nM8R3P2+t0ki0g3IOHDzxhulrtiBLWImbg1LKfh1HU3l3SR3buEYXgFw3/VPb7fLyehziJgo4gsJ9CB3zl2CmfZ4eembkF2kMjnbLY+T+e5vmzfJfRZuk/GsJPaCdyFIajm+qzedwAXJeFWv0qkcNZ4HM3mzph+Euscftcvt7S/N41nKLNbw4pg6W7LAj2G+7yXSPmsiUuYzZ2xKYj38zZrVWZ/VanWtpw3ADsLEflAVU8UkTIChYcAqqq2Qh9lkYblzGlstrN09pqA96kUXE5cm3UKqnqi73frROrJosPPbmd0CmbnesWzRqXlnk7rz8lGRNphxAo8DkMDukZVPzZdPxe4G8MnVj1GOPsP0t3OuM8SuqQPpz1LyVrATrczvSwf0+YAi1R1CDAS+DTo+tvASFUdBVyDIdzSTgJBKBz2FLk4ErPQSqUAy+ZdQhFpA5yMEa4eVa1V1cPmPKpabjL2bEkTLWIncJbQxSU8yegjfv/q6dG0nLtLmARL9/7AfuAJERmJ4UXhBlWtMNcjIucDvwM6A9OT0vgYkUiakojoxo3LQ0JwV1XVUl5eRceOrS1Dd4cL0934OzRfYB2JpZuvmamsrHFctOqqqloKC/MdtQZXW1tPTo6Ql5cbV/nD+4ywCe0629o48gkaCXnv/8j8n12k9Lo6D7W19bRqVRgxf/g6Aq8Fvg6fR0Q49tgJqGpcX7CI5HeaSO25S8PneW0q7H6H9sFak6mOccBS4ARVXSYic4AjqnpnmPwnA7NUdVoc7b0AuBdD6AnJ3iX0erXhg/X/5OfnkpOTQ2FhQcgX4v8d/guNnmb8hvCdLny+wNeBi+0rV37O2LEDov3JGcXatdsYOrQn+fkJWKCkmW3b9tGqVSFFRfYETqIED1yNr6MPdv7XpaWVHDpUQe/eRWHLWQ2+ga81pGz4gdz47fV6k/pZxMlOYKeq+oPv/IcIrmR8nh0GiEiRqsbq4/33wAxVDV4js0XUp6BHjw4haWVlVdTW1tO5c9t47ukSAw5SrJqM4MErHqqqasnPz3WcBp4MVPUrEflSRAar6iYMY9GN5jwiMhD4XFVVRMZgBKQ4GMft9sYrrCDLjua4uEQiyweHnwDPiEgBsBW4WkS+D6CqjwAzgW+LSB1QBVwcp8eFFSLyPEZEHvNZwhfsFM4qw1GXdOKspz+bdwkBVHUNhusXM4+Yrt+LsfaUKG2ASuB08+2BVAos/z1cXEIxItA0dStiI9t9uqcLVb06kfJZ49PdJZ3E1jmsgps2DdktsUQkV0RWi8irKbxHTxF50RfwYq+IzBeRnnbLx2np7q5huSQHv61VU+NErTAF3ECohXuyeQLDz1Z3oAewwJdmi7g1LJdM0gycSyYIK4PsHoB9Ws50Un/kppOqPqGq9b6fJ4FOdgu7R3MSYOnsNY4UWh6Ph5dfXkx5eUX0zFlCFkR+vlJEVph+rgu6/gBwK8bh5lRyQEQu900/c0XkcmIwj4jTrMFdwwLDX9PS2Wsc5atp7doNPPbky3Qachpv3vkg3770DCZMGJ28G+DU6ZVzO/T+8hY8umxI+AxHNgNlYX26i8g3gH2qulJEJqekkY1cAzwE3I/xoX/kS7OFe5YwQW7jNlLb2ZPzYZeXV/DQw09zxNuZ4pl3kpOby8Dhk3j5v3P58KNV/PhH3yY/Pz8p9zJwVidxppBNGicA54jI2UAh0EZE/qWqlyf7Rqq6A7Dlv92KrPGHlQpm3nQGB2c3dSuis2jRu7z+9kqOPeVq+nbs2pAuIgw/+WL27/6cn95yLzf8+FsMHNivCVvatBzlU8KwqOrtwO0APg3r5mQLKxG5VVV/LyIPYjHCq+r1dupxp4RHMXv27GXOg0/TutfxHH/B7WHzdeo+gA4z7+Rv/3yc4QNWcsUVMxN8eJ3XOdwBOOX4dx9XJFKJaziaIB0PZd7n4PF4ePrp+ZRs/poR026ksHnLqGVy8/IZe+b3+HLzSn5++738+q4bad68MK77O3V65cQ2JxtVfRcjJH2y613ge1mpqvPM10TEdjQe13D0KKOk5FNuuvleDhWMZMKM620JKzO9Bo2l/4nXcdevH8Dj8UQvcJSQ7Udz0oiVqh9e/Q/CPfyc4dgd9SsqKnn4L09zqK5Dw6J6vLTr2JUeoy/mnnv/yi9u/1GWrO1k79EcEekFPI0RkdmL4dxvTlCeyRiBVL/wJb2gqrZXcEXkLOBsoIeI/Nl0qQ2Gj3hbuIefjwLeeut9Xnl9KcdOvpreRd2TUmeX3oOprTzEP/7xLNddd1kcNTjr6c/y/lwP/ExVV4lIa2CliLypqhuD8r2vqt+I8x67MdavzsHwaOqnDLjJbiXuGpYjsP6s9+7dz5wHn6Kw63gmXXhH0u/aa0gxm5Yf4KWXFnPeefZtzZy7huXARicBVd0D7PG9LhORTzGOzQQLrETusRZYKyL/VtW6eOtJYEoYb0mXRPF6vfzrXy+w+rP9jJx6PYUtWqXsXoMnfIOlbz9Fp06fcMIJ41N2n6bG7hJH8MmGzDlalBxEpC8wGlhmcfl4EVmLoS3drKob4rhFXxH5HTAUw+YLAFXtb6dw3FNCV8NqGjZu/B9/e3QevcdcwMQZI9Nyz5FTvs1zL/yeAQP60LVrZxslnNk3HKtgVXSGVd8Jf71sDrApUhAKAESkFTAfI+bgkaBaVgF9VLXcZ2D6EjAojtY+AdyFYel+KnA1MawfJOCtId6SLvFQVVXNfX/8O0+/uIoJF9xJjwHpEVZgDFBjz/4xv7/v0aN2syULdgmfUtVxpp9gYZWPIayesfL+qapHVLXc9/o1IF9EiuJoR3NVfRsjAM52Vf0VMMVuYddbgwN4772PufWOOTQfcC6jpnyb3Nz0B6QobN6S7iPO5amn5kXP7Ejs7RKaw44dLdNBMR7ox4BPVfVPYfJ09eVDRCZgyI54fLpXi0gOsFlEfuwLHWZHbQcS0LBcUs/+/Qf5y1/+wUcbq5l04S/p2KVXk7an56DRrP+inM2bt0bN67QxLVYN62gRVj5OAK4ApojIGt/P2SLyfb9fd+BCYL1vDevPwCVx+nS/EWgBXA+MBS4Hvm23sHNiR2URXq+XZ599iRXrd3P6Od+iru1QvBkiAUZPu5Y5D9/N/X/4edjD0k6dXmXIR5x2VPUDonxhqvoQhpeFROmrqp8A5RjrV35Ld6tF/hCOSg1r/v2LHemnCmDTpi389JZ72eM5honn/JQWLVO3AxgPuXn5HHPit3no4aebuilJ5Whdm8tAmsbSPVNxqqD67LPNzP3PIirpyPjzfkFuXjJdvSSXzj0Gsm7zJ3zyyRrGjx9lkcOZD7+7Nps6kmXpflRqWP71hUwXXl6vlzfeeI9bb7+Pp1/ZQN8Tf8SoqVdltLDyM/yUS3j8n6+GjVzstGc/VRqWU7R9EXncFxhifYQ8k33rWxtE5L0Yb+G3dK/GsHT3/7wC2F4QjKphHT5suNENDrVdX+/hwIEjYUJ0hw/nHRwu3G5I8dAw5Mrna79siPbbbE0LEKgSZfELy2nVtz29h3Rn48YvG/6WqqpaNmzYEWCSERiyvOFV2DzB+QLfB+QKa/rh8Xh444232L5jJ8NGjOO0M8/3je7bIMgGOFdr6FC3JWkSIBmPpb8lnuoqVq82jpaZoy5XV9dy+HAleXk5QSUC/wyzRtP4UgLeh8tjVY+/DXauBacfOVJJfn4ulZU1FvkDy5nTrPKZ6z7liuMBOHiwLEz5jJHsT2KsUVnO9UWkHfAX4ExV3SEitnf2IMDS/RlVta1RBRNVYB04cASrL8Lj8VJRUW1xrTGP/3VOjvheW3eYxi8tfIdqTG98v/HgZ7T9R2dQGKgDEa8wSAex+cXNHP7BXvqc0SmgTEVFNQMGdG241+uPvsvZ350c8LAF9h8huD8Fh0OPlt/8twCsWbOBp555hX7jL2Tq1GFA5L3h9nWfU5rXG69kltZ1cO8OJkw4jjFjAg2UVWHz5t106dKONm2ahwhta+FvJfjNg5ZVHutBxc5AZ5WnpqaOZs3yad26ueXgGS6t8VpwvvCDdmiepkdVl/is3MNxGcaB5x2+/PtiqV9E5qrqN4HVImLlwG+EnXqiCqyBA7tZpn/9dTl9+sQkZJOGX8XuOLsHwxkOwEVjDjRcn7UKPNX1tGjRLKBcTk4OzZrlIyINdbz6l3fSskVdWlrGQw8/RVVuT8ObQk4ss/EM6dUmtq1+jV/ceEGIhuAfgHJycmL8G5uWsrJKWrYspGPH1k3dlNip6AKrfhD+etl/sGPpHoVjMIxF3wVaA3NUNZadlxt8v+M9PA1k2KK7nRh18+9fzNLZa5jOdIYzvEFQrR4zEYDRq8LvjorACw+8wdLZaxrSimdZLRonlwUL3uStJesZdurVtGnfNELAInCDAAAgAElEQVQ+mXi9Xgq8pbRv3y5snsyZ6dgjUzSdFBI2CIVN8jDspqYCzYGPRWSpqv7PTmHfAWtUdbuIdAUmYIzEn6jqV3YbkRFDoHlhMtICpf9aOGFlxiyUGgl8ivzCKtWLoh6PlzYtctj0wTNsKfkQT30sh9Uz78nf/tlyTjlpTFM3I6lkQZivRNkJLFLVClU9ACwBYj4fJiLfAZYDF2AYoy4VkdRGzUkGdgST1W6foV31CMg/etUyVo+Z2KBd+aeJ4TS2YK0q1VPC8847g/POO4Pa2lo++mgFHyx5kMPlHgo7DqTv8Mm0atM+pExFWSn7dn1O8+7NDAU8gziw5QOmXPnDsNedadPkxDanlZeBh0QkDygAJmIcYI6VW4DRqnoQQEQ6YoT6etxO4bQKrFDBMx2AhSxsSDcLk2ChZp4KAsxbZZy9vGjMgZCp4HCGUzIb5hMotM6/4bQmW1spKChg8uRJTJ48CVVl06YtvPnWPLbuLaWqPp+8XKUwT8nPgw7tWzOyX0/mPreI48//KW079WmSNgdTXVVB5/YF5OVl1GpCwmS7hiUizwKTgSIR2YnhUSEfQFUfUdVPRWQRsA7DK+mjqhrWBCICOzGc9vkpA74MkzeElPc6K03KLKyAgNf4nK4uZGGD8Fo6ew3Fs0YZ7y2css5bVRSw6H7RmAPMW1XEcIbzLrsa0jOpP4oIQ4YMYsgQw0NHXV2d5VGXLl1689zcp+gxaiZd+w5LdzND2Lzidb597rSmboZLklHVS23k+QPwhwRvtQtYJiIvY6i15wLLReSnvntYHr72k1JVI5ywMuPXliIRPIUroaTht/+fX9vy4xdggeHkM9e1c7hzec2aNePXd92EZ8+7bF373zS3KpT60q0MHjwwar5MGhzsoJq9Pt3TzOcYvrT8T+LLGN5OW2Nj8SNlAivcGpVZ+PiF1fCgf2BoXWbh5nfrcXCWoTH5hZZ5OukXWvNWFYUIsPn3L+bgrkOJ/llNQk5ODj/76Xfp22Y/65c812Tt2L9nO0MGRvcZ79TDz85sc3KwY+meDFT115F+opXPiF1CM35BBI3CLXhR3C+kFrKQ6UxnIQsbypkFlbkuPy8/9FbS25wuLrvsPKZN7MYnCx9ukvt/vuw/XHDBmU1y71TjVD/0SeRJIGVfrog84Pu9QEReCf6xW0/K1rDiOc9nFjALWcgfbvq5Zb1LZ68JEFZW5c10nO3bVbzbdlMyDEN7rq6u4b0ly+nYK/2+1b8oWcIpxYNo27ZN2u+dHjJ0rSBN2LB0T5R/+n7fl0glCQksY94feVgya0fBwquEkoYpYLCwimTQWTxrVKMQ8r1fOHth4OK9D3/acIZToiWO1fq3bdvBffc/zdAp19Ghc8+03ru89CAVOz/i/Ot+arNE9PWg0tIySko2smbtZxz4ugyPB5o3E4YNHciYMcPp2bNb2nftjvJdwkQt3RNCVf2hvVYAVarqBRCRXKBZ2IJBxC2wEvlyzbt9VoKqmFG2baP8gs1fp1/zMjOc4ZRQQl3fGmacOTXudkOojVg6WLz4Xd5dupnimb8kL78gbfcFY1Ba88Zf+e1dEY5+xMATT8xl0xf78ea2onXXoXQf9E2GtG4LgMdTz+btm1j2r4+pLd9DiwJo16aQ0aOOZcyYEbRtmzqDNGfajhm0qChnSIQTHpvLjlCWuKV7sngbmIbhwA8Mq/k3gEl2CqdcwzIToG0ROlW0I6jm37+4QbtayMKGMvPvX2wsyFsIrXu4h+JZoxjYezCxqv5WU9p0Cau6ujqefOoZ2vQczcRzbMeaTCobP5zPZRdOiUlYhHv2V65cx5Z9+Yw6y/pvyc3No2f/YfTs32i+UV1Vwcefl/DM83/g4Tl30KyZ7cE4Jtw1rLRR6A9mAeCLwtPCbuEENKzEzl/5hYz/tV3Mlu7BR2vMu4p+oXVw1q4GobZu3bb4G5xmdu/+inv+8BijJ06hVe/j8TRBG3ZuXknn5l8zadI5cZQOfPpramp48l8LOf7CO2OqpbB5S/oPK6ZF6/Y899wrXHnlRXG0pRHzABTc747yKWGmUCEiY1R1FYCIjAWq7BZOcEqYmBodq6Yy//7FAdbxxYSuc5nXvg6yK+gesdlhmb1C+M0p0qFdvffex8xfsJRx59xGl9w9BAeISwc7Ni0n5+ByfvzT74bNY+ewup+H//JPBp90VdynDLr2HszSF1/k0tpaCgoap8WxasDmgdKMk6eEycDK0l1VH0vBrW4E5onIbt/7bsDFdgsnOCVMpHRsmDtZ8KK8WbPyW8VbdVpjAI2t0f7pp7/eVKGqfPbZFl5Z8BblOT2YNNO3Q1oH6d7B2vHpx+QfWcONP/1uRK0jkmAwF1u7dgMHqtsyoltix4v6jDmX559fwBVXzGzw2gGhhsXRBOlRFvEmKdixdE/SfT4RkSHAYAw1/LNYQtcnpGE1xaDkn+L5CT6TGE5Y+bHb5pC1sij1xoqqsm3bDpYsWc7WbV9RXgMtigbR67jL6NO+U2O+pN3RHl+s/4CW1Ru4/sZrk1Kf1+vlsadeoXhmbFNBK7r1OZalL7xEwX2vkZebm9QBJNvPEqYL33rVTzGiSH9XRAaJyGBVfdVO+YTWsOw+TsEqeKpGuFRpQMkUVitXrmXBwveoqBEK2vam+zEnc+zpTRtv0M/Wkvdo59nMD39ydUL1mKdXr7zyBj1HnJ20A+d9x5xLnWctefuTHU0ou4/miMiZwBwgF+Ng8z0WeSYDD2Acij6gqqeISCGGq5lmGPLkP6p6V4RbPYHhy/143/udwDwgtQLL7npQsnxNWQmMW9rf2/A6nKGpmVhH0IOz/Gtgx8VULhhV5d13P2Lhog9p0WUkg6fcGEP05vQ8RZ+veYei3O187we2Y1pGxev18u4H6yie2bho7/V4OHL4ALXVldTWVFFTXUl9TRVaX42nvgqtq6a+rpL6ulqOPeGbtGwd6CSwa9+hLH/pFf50b6BNWHA/i2V9DbLCgV9YfLZQDwOnYQiQT0TkFVXdaMoTzqd7DTDFt9uXD3wgIq+r6tIwtxugqheLyKUAqlolMTyYCe4Sxv4tJ1O7Kp41qsHq3a52ZbfJyWinx+Ph1Vff4r0P19Gh3yTGnntHRk47tqx6k26Fe7j22m8lrU4ReOmlxXQfflZA+uq3Hmdgtxy6tGtHizaFtOhaSPPmhbRsWUTz5oW0aNGCFi0K8Xi8/N899zPwhGvoGLT21axND/bs2Uu3bl0a0qwW0+16sAUYeuawjPxu0sQEYIuqbgUQkecwvChsNOWx9Onui/7sN1PI9/1EespqRaS5P4+IDMAQerZoEsPRZBLuvKEV8Sy6x0NNTQ1z5y5gVcl2ug09jQnnz0jp/cqPHKKgoJCCwuYxl/3fitfp1fogV191SVLb5PF4WPLxeoovOLchTVXJ93zN975nz6bs9/fcym9++xCVZVPpdcy4hvQOfUbxySdrOeec06PWEU5oBQu3Lz/bQ58+nULyZQk9CPRJtRPDQZ+ZsD7dfRraSmAg8LCqRorifBewCOglIs8AJwBX2W1owoajkUi16+FM2+2pra3liSfnsWnrAXqNnsHE8y9Ly32/WP8+S17+G6edPpW6eqG6Huo9RrSiHBFyfOGoBCUnx5hk+tOOPaYnl1zyzaS2RxVefvkNegw/OyD9yKF99O/TNUypUPLz8/nVXTfy10f+yaZP9jF4vFFf977Hsm7JmwECy8oWz2ynF+5sa8Pur+DYNaxBlXnMDfJOYuYK8vk48tEcq788+OEO69NdVT3AKN+08UUROc7KuZ9v6vcZhnvkYt99b/C5XLZFWncJm17ApG5nc9my1Tzz/GIGHf8tJpw7IDU3CcPwSedQ2KwZnQt2873rrKd1sa7pJILH42HZys+YdFGgICz9eh/HdY8tCIeI8IPvX8Gtt/6GzQWtGTTyJHJz86jyTSKsBsWGXePZC0OMi82Emqo4VGLZI9LRnJ2AeeenJ0bg0+A8B1S1AsP40+/TvSEIhaoe9mlgZwIhAktVVUReUtWxYPJaEANxb91EWsMKF+22qSPgpmJKWFpaxt2/eZBX39/J8RfeSaee6RVWfgaNPYNDOYO4749/D/lezAE+0vEdvPXWO/QZGeqppLJ0H126hNcErFi1qoSf3XovbfufysARJzakV3vyqa4OXPqwMnGJhPl6j0FdHKthJYFPgEEi0k9ECoBLMCIym3kZOElE8nymCROBT0Wkk0+zwrc2NQ1DiwrHUhGJ291IytewUm1wGSvJ1LAWLHiTN5eUcNyUa2nTLkXrHzE8RH2Hncjuz1sz++4/M7TlMU3iu76+vp6Nn33OGVddHXKcqLbiAJ072xPoW7du5x+PzSOv/WDGnXsHObm5AdfbdBvGnF88zsBejYFci2eN4qDPJbbVKQgr/Frn+vU7bOU/GlHVehH5MbAYw6zhcVXdICLf910P69NdREYAT/nWsXKAuVFsqk4Fvi8i24AKjB6uSQukGgkrDSt4BPfv5GWC0ErWRsGuXXuY8+A/advvJCZdcFtS6kwW3QeMZP/GZcz76EUumHoO+XmBrpdTOS3csuULHn7kWU48dTpWkra67CCtWkU+57p//0Ee+fuzlHvbc9zpPyW/oNAyX89BY1j+3isBAssOka3zs1fFUtXXgNeC0h4Jeh/i011V1wGjY7jVWdGzhCeqwKqursMqxHZ9vYeKihq8Xg0I453XqgBBGvurwEn3FTeo2wcOHAmpy1+2MYy4dTjx8OmRQpM3Xisrq2LHjv3k5eWGzWvOH5hu/M0LFrzGnn2lTDn7IgoLC6H+i8APTANfhH8E7Kl6+VpJW92B2lS1cg7vp2PxcDoO7MrLr7/KZRddarQT6NK3iHXrtoeUsfOcWj3M/qS6ujrmzXuJqjph6tmX0SK3lrr6L1HJwfwJjBs/kZ/dci9nn3kqxx03LKCOyspK5r/wCkcq6jnptAto2aoNsB/qaahDgdqaGtav+og92z9j3GljadHG51Aw6PsbNmkgIGz4cLM/keEnD+bLLw80/C3GL0EEqqpq2bevlPz83IY0899tzuv/HS2tsWz0PEc7PgPT72PsJJYAj6lqfaz1RBVYW7bssfzw/V9wQUFeQPqwU4cE5N30ydaGznRs8QAqKmp81wLrDE0DEHJyhBWLfet3qow/a0SYThT42n+t8TWUl1fTrl0rCgvzLa8H1hVYx2ef/Y8nnn6JQRPO5+RThlGH75hfY+4In6JEEU/hy7ap30F5ble8Ys8PVt72TQAUksvUq37NosUPM+uOH0aM0gzw6t+MABfTr5scNo8GCWOAZctW8fz8tzjupMvo1KMfVUBB3Q7KczvjlUbtToCCTl0466q7Wf3+fLZ8/grXffdyROD5519mw+Y9DD/pYjp26g6o6fi+Ul9fx+Z1H7Nv2yraty5gyqmTGP6ds/nghZXUldX4Pr7G7/L4c4wB/+MFpoAnQkNfbfxbQgdij8drMQhGG0CD67Gu21yH+VqW8BTGI/M+hpY1lMbw9baJKrCOO663ZfqmTbvo3r0DrVtHtv3p379LxOvRCJ5idukS+cGLdAxo797DtG5dSMuW1tMMK6qra3j4L0+zv6o1k775K3Jz84h5WEgEyUElL+Dhj4TWewGoK55BK2DcObfy69/ex523fzfA0DKYC35yWkzNOny4lDl/fhJvy4GcdMmvEZGGNSuVXDzSzFrI5sLwyVdwYPcX3HX3g+Tk5tFv3ExOuPBYo91B2Wurq/hw7my+ffk5HP+jmwPW5erKrO0Ni4raWG4uROo7Bw6U0bVr+wChlk3YOZrjyzceWApcrKr/iaHsUFUd7sv/GEb055jJ6MPPZtuaSGtg5nxWwVnBL7hiU7/ff38pc194lyGnXMWoLtaCO5PIX7qAuuJAI9XmLdtQfMEvuPve33Pz9ZfSv39iHhNUlRdfXMR7SzcxfMq1tGrb0SpX1HqKuvej6Ju/ippvw/vP8ovbvkefPvbPW8bndDF7zxLaOZpjyncvNHrftFsW01jkW+SPq60JDifJl1hWnS2csDLn7Ti7R0AI+2BbnPn3L6ZVn/b07Gn1gAXy9deHeODPT5HTbiiTLrrTMesMwcLKT0GzQk6YeTsP/uNhTpk4gPPOi2/hffv2L3nw4WfoeMxUjj//3OgFEqS6spyWOYcshVUs5hl2NhqcG5osKdg5mgPwE2A+MD6OsiNFxO/aTYDmvvf+XUJb0U0ySsMK1wmDO5x1CPseWDGd6Y0GhApv/esjzvv+FMu8qsoLL7zO+8s3M3zqdSEHb51Mbl4+E8+5kQ1r/8vKu+7n57dcR6tWLW2Vra+v5x+PPsvnu2sZOf3WsDt3gUjMLrSD2fD+s9z43QtD0q2+fzPx7Eg7OZBqFVVhI0YBlBtH/SJZukc9miMiPYDzgSkECiw7x3pQ1dzgtHhIi3sZO8TrOz1Ys7LCL7SmPXpq2Dzl5RX8+u6H6HjMNIrPOy96g9NK8j7nASNPpbTXUH7+yz/znStnMHp0ZE8Uq1aV8MQ/FzDo+MsYN/IYW/fwejx8tHAOLQvqqNMCcpp3olP/CXTrbd82rLK8lPbNyunRo1tAutURnOCgI2YLd5cGIlm62zma8wDwc1X1BA1CdsomjYzxOBrOdW0kjE4bvuOGYt3gysoq7rjzfoad/pPUGYBmEG07dOGEi+5k7mv/Yumy1Xzvum+FCJIjR8p48KGnqMrtyfEXzrIUNFZrZh5PPfOfe5Trf3QZI0YMBWDPnr18+OEnbPzvQo7UNGP06d8jNy/yJsKn7/+bm3/UeLTH3DeshJU/MpI/zW9AGlw23CCY5Q787BzNGQc85/uMioCzRaTeZtmkkeCUMLmC1G5gCn++4FF04WxjoT1kxPUd07CKjFJTU8Mdd/6JY6f+0JHCykpo2CEnJ4eRU77N7q0l/OyW33Hrzd9p2EVcuPAtFv93LcNOvYa2HUJ3FvOXLgBC18zq62pZ+uLvueDcsxgypFEb69atCxde+A0APv98G3/88/8x4dybKWxuPSWtq6uhdUElXboY30ckrcof19Iu4fuXc6eESaDhaA6wC+NoTsDJfVXt538tIk8Cr6rqSyKSF61sMmmyqDnhsGuJbZUv3Il8v2O///1vN127tm9Ir62t5Y47/8Sgk79Lu472vQhkEvEIKzPd+w+nY/cB3PPAXyge1YvVazfTru+JTJp5e0z3rKut5uMXf89tP72cw4fDx/gZMKAvv/7lddx97xwmzfxFwDW/IFz36Wpmfs8wiA723e7fCbYSVGYtK1ayyyQqEDtHc2Itm6q2JuRxNFMJJ/TM62719fXcOetP9Cu+Ku2RlDONZoUtOP68m/li43KGnTGDZoX2wsT5tbua6kqWv/R7fvnza+jevSuHD2+LWK5Tp45MO3kEGzZ8SL9hJzTU5eeLXdv54r09bFvyFdC4iB5JWEFjUN5wEZXMWHmvyOIpoa2jOab0q6KVTRVp9ziaCXg8Hu761QP0HP+tEG+W2Uy/oRPCXjMLFD91xTOoqixjxYL7+NUvv0+nTmaTEaWsrJzS0jIOHy6ltLSUwYMHUlRk5Jkx4zT+e8s99B4yMcBd9OGyUtq0ahMiPMw2dpEItr+LRKDQcmZfzjYyyqwh9Qgej5fZd/+ZzsNn0ql707iCcRJ+QWXlGA9gwT/uoE+P9jz01+fweIV6L5QdOkh+i3bkF7amoHlb8pq3IaegNS+8No/+PZpz7TWXUFjYjGuvOp9/L5jLiMmXUVc8g/ylC1i+YSUTho0NaYdZWPk1Kb+m5X8fHFHJDv7lg3ZDOvPCA29kgM+2piGatbqIfAvwB00oB36gqmt917YBZYAHqFfVcaSIjDMcTS3KH//0NzoOnUGX3oNtlzJrF4muGcVO001TrLQqs7DKX7qAb008LeQz6VC3hUN5fVEJ6l6jTuTAnm38/JcPcvLxQ7nggrNo/tKbVFeWU9jCiILzdekhamqt/FyFmq6Y16vsTAPtkE5Hh5mCTWv1L4BTVPWQiJwF/J1Ae6tTY/EcGi9Zo2GpKo8++jTdh55Gl77xR8EJt0N2NOLXevzTMX/06+iE7xhF3fpSNPMXbNn0CTf+7HcM6NuZHZs+4ZjRho3chdPO5c2l/6VFYXNOGjOJHMmheNYoFs5eGHZKGEsQEjMB08w/Nh4Bi9cm0MFEtVZX1Y9M+ZdimC+knYQ8jjpJw7r/gcfoMXAcvQYmFrIrWzk4a1eIRbkfK00sGn0Gj2fizFkc9PTEW1cNGAIyLzePs044jd5de/Hc4vkcLP26ocxC0z9zWqwsnb3GZ3AcKACnMz3s35gub60p5EoRWWH6uc50zcpaPZI19rXA66b3CrwhIiuD6k06KY9LmAk8+NCT1LUZxeBjetuPJ+TSgF97Sbb1eE5ODscWN8YsNAu+fj360KNzNxZ/9DaH3yvn8jsvbFiI94d2C8aO8bHdxXtzfj9+7SvTtK1DHIoouA9yEBK3dDcyipyKIbBONCWfoKq7fbEK3xSRz1R1ia3Gx0hKfLpnEn/7+zOUFwyhz9BJZLIpRqbin/qG0zz8hGpZiX/WBfkFzDjlLCZeOZrn33iB0nLj7Gyw4DS/T1SY+IXZ0tlrGoSb/1+0z8DB2LJW97lDfhQ4V1UP+tNVdbfv9z7gRYwpZkpwfFzCSDzxxFz2e3oxcNTJCdWTDetVkagrnsHY12aE+KqKZyoYD8f0GUjvbr1Y9OGb9OzSg3FDGz3yWsWlDGdA7M8faT0Mwp9P9Z9JDb7fUUBUS3cR6Q28AFyhqv8zpbcEclS1zPf6dGB2qhqaUJSCTNawnnnmRXaWFzFw1LSgK5nb5nBk/tCQOOGEn1+rKSxoxnmnfoM24ycz780XOe6WQRGD6Npdbxru+xctj59YppNOweeq2G+t/ilGIIkNIvJ9v7U7MAvoCPxFRNaIyApfeheM8PRrMZzyLVTVRalqa8YdzUkG8+a9ypb9LRg8MTjUVDY8+unHrIHmL11AzqBO0MF+cAi7mlpd8Qz6AT0O7OHV9xcx7thRDAgKQhHvwnhzmlNlcswMsZ9TdDLRLN1V9TvAdyzKbcWIT5gWEpwSZp7EeumlxZRs9zL0BOtpXDJEVvADls4po9Vh53BW6KnG8kxh8Qz4eqntOoJt3MzvzZ5mzfcqOPVizp/8Tda+/SQDehvnFoMPNUcSXOZpod+WqznNs0pAOZWjSsN67bW3WbGpguNOvsTyupr+N2NHAEXSAuL1mGAHc2v9bbCjkaSyTdHuyzGdbEX5iSSsoiEijJp2NRtXvsGXcx7nhhuvjrquGo8Lo2AWsrDhML1L+knQrCFzJNabby7h/dUHGDnlipBr/oc3d9dmvBW15B+pbniYrR4SqzT/jpF/+zh4pyrcwxbpPsF5wqMNBpx+zPe38nnfZMIqBoKnkmbCaVfBDBp7Onu+KOH2O/7AXXdeT/Pm1t5Qoy2SB08HzZgt6o9Wx4A2juYMAZ4AxgB3qOp9pmvtMHYPj8MYY69R1Y9T0c6M8TgajUjuk99772Pe+ngHo0+7JuR6gFbSrU2DzyOrh6vjbGNnKPiMGsBtBK6HmR3EhRMk5vuE2xIvnjUqqja04XuPcMrOUwLaUDK7JMD2xu4DniqSvWMYy9/Srd9wmrcu4tbb7+OMaRPxepXcforH48Hj8eL1Kv/5z0L69etJZXUlLXzeKBYSebcQgjxAzBp1tO0OAraP5nwNXA9YueOdAyxS1Qt9oe7tufuIA8cfzfnooxUseOdTxp35/YD0WB6gcG5Lwq1pBB9RMXd6syALrt9PsKABKJ5l/XDm7N/JKZwS0pZgQ8HgM34Q/WFPxbSxrngGWve/6BktyiUi9NoVdWP8eb/g051bkJxcRHLIyclBcnLIyc/FU19Hycc72LP6Y6qqq+Ak2PX+XtaxjmlMY6w39MB1sGZ1NAorH3aO5uwD9olIgIQXkTbAycBVvny1QG2qGpr0UPXpZOvObby6ZBGn3/IQ0PgAhuv4Vu3129xYCSezVwCzR4Bw0zy/IAue6vgXecEnWGYTVoCZyy2dvYbpN0xnkAziojEHmLeqqKFd0aYmdgSR+e9IVHAFlo+9XyRDQytoVkivAeGPXnXrPYgRdZWNCafCktnL8eLl/dXvk0OOZT9I1sHqVNG6R0uKrw3fvg+eXsLubYkFoYhAf2A/8ISIjARWAjeoaoXtPyAGHGs4un3Pl5Rs2cA5p5xNvUhMC9J+InkBCD6vFm6K4q/D36mDr/vfF88KvX+44xTmA8ftMTykmoVVuAcoVqFj/qyaapG+qSmggIlM5MndTzKJSSHXYzlYnYnHdkwk5WiOBXkY61o/UdVlIjIHuA24M4422rpZRLxef+huf3ht47XH46Gurp6amjpTGO7APMGhuo3fkcN9B+f1pxe0b86Es0bwyaIStu3czpov1nPRBRfiHTCKfM8B8opM/sFNocuRxm8jv1UBuc3yyC3MB4GBFw00LgqUSAkIbJbN9DypK1Ok0Tre231gw+vC+l34v8vcnm0Zdu0wtsk2pk6YjLdrX5rXmwcqPwrjRjY0bcQPvdSMqWAap+Lt3JucfTuMbLs/AIH9h5TpPzE0sENDDuG53MORiiPsZCedRnYwynWyCCxav93i3uHJ6dOe/WuMw8WdRnUw7m/C29k6eGxgTw7t67laSxvProZrkfMb73O6twl5RFSVZp6vguoILR96zXxfCcib2755Y36Fk/90PCufX8m0U6bx1hdvMahoENt0G3hhmXcZJ11uKBonXzSe0tIKQBBpHLBFhLf++aHxdzfLo7KyBhHxrZX682JZzvy+iUkkkMROYKeqLvO9/w+GwEoJUQXWmjXbTF9A4wddW+vB6/VSWVkbkB74OvhLakxr/B1aNjQvjJo6lOrqWjZ8tYEDVcoFP7kLzcnxdVdpEG4An7+0A3yCceB5fQyBCuTWK956L1uil3IAAA2eSURBVJ4aI9h89+OLWPH8ClAYfHE/6gdPoIgBGL4DxGKICepc3cfQ+nzjZUW4PKZ0f33Nzz6NMvOlrkXklSxh3d8+42ROpqe2YAlLGPH9IbRo3prtk/rhkR4UcmxguQTIW/cu6/62iZPxCWa/siVw+Lp9ANR3LbL4KyDy4KsUeMuoknYgEiF/oJjJLW9c9tj0/FYEYfAl/eB/JQEDkGfQ2IDy1vWrbxxS0zXfIDhwRMPr3B2fgkCr3q3o2bsn77Z8l439NtKioAXkQFFOEfltCxHgwAHjHGPwILxt/S6ad2nt+0OED95YR49BXbEasM0DsvlaBhD1aE44VPUrEflSRAar6iZgKqFBVJNGVIE1Zoy1xfKBA0coL6+mb9/OSW+UFTt27OL+OU/SZ/w3Gdt3mE+oNJJ/sFFk7HxnZ8PrXic3RsLJLcijtvdIqnPbN0yH+p3RjbriGaRkwh0DUlXHydsatbphswZQPfJMmtVvxyOF1Oc0T+79qusp226IP/O6TQkleKqMU4P1OdabPdHs1lRyqctpFRqiKAL5R4xvNGCD4hPjV8CUPLe9VfG4yN+7HIAt87YweOJgztp4Fos2LuICLmi4Z+WuUgAGfrObZR3rX1vf0G5/G4edN94yb6ZiJwiFiHQFVgBtAK+I3AgMVdUjGBGhn/HtEG4Frk5VWzM+CIWq8vzzC/hk/VeMPfd2m1GHDYLXnQo8exuuZdJ6TeOaVY+4ndHFQ6zn4lJtUR/O9CM4zmCy8G/QFM8aRamUspSlbGQjU5jS0B4730WDg0N2ZfIaVkRsHM35ijBO+1R1DUbcwpST0e5lDh78mtt/8Qe2V/Zk/PQfRxRW5geneNaoCHY8GaOGB1A8a5ThkzzIPCFnv9WaWOLUFc+w5fwuf+mChh+/47uOs3s0CBfrTY7YPmN/HVbCYSELA8xIkm3v5e8frbq1oHjWKIYylNd4LUCYhxNC8+9f3OAI8OAs5worJ5GgwEpmUwJZvPg9Zt/7NIOnXk+fY4vjqiOTtKhI1BXPCGirWdOwXFxPEuG0h+D2AAEeOocznNu4rcHQ1pIYpoPB9/ILKbMA9/unSgV1xTMa+nIzmtGa1uxhDxBZWPnbM/nQcY4XViJypohsEpEtIhKyaC4Gf/ZdXyciY0zXbhKRDSKyXkSeFRH706AYiTolvOOuh8nPU9q2bs6550yjf39/WKzUHH4uL6/gT/c/Bu2O4/jzb42pbHS7IknTRNYZdJzdI8Cq32xjFpwPGte6zDZhxrSpMW88g4RZawoWoo1rWo3Ts2haVrg22DHdmM50aqjheZ5n+qyplnnMpy4mH3K+y22blu5nAYN8PxOBvwITRaQHhgX8UFWtEpG5GIv2T6airVEF1vAzbwKgquIIj859DU/ZC4w6rh8nn3xCUjWsysoqnn32ZdZv2sOQU66yDJFul0zRrGKJtmP10NYVz4jZVCEWDs7axbv4hdQY309gm8zC6qIxRlCUeauKGu3BLI4iyeDO5G+KbG1vPhtpNX03TxMPsouxxXc2OBCM5Vym3YPiDO7csI7YjGYc5CDb93wZ1rZq5k1nZIxWVVZZzvL1K8JeP1JxJFoVUS3dfe+fVmMdaKmItBMR/05EHtBcROowjuXYNYmIGduL7s1btmHEKYYXhD3bN/GHOf+GuiNcesl0JkwYHbc9yeHDpTz19Hy276mg75jzmHh+6mIFhvPWkAh2jsHEKkDTJXDt3scvrMzGq8ExAMOdMAjnDieSsIqlbcH1xHuIOpiBDGT5+hX06RY6Jc8UQeWjfv+hA9/5+wtPdoiQ5whwsYh815QWq6W7ZaAKVV0hIvcBO4Aq4A1VfSOeP8QOce0SduszmN69fkB+7QEWf7KZuS8+QL9eHbhw5ll07WrPzGH37q946ukXOFiey8Dii5gwvms8TWlSogmrWB6QTNEKg/FPAyNpVRB4yPvE+05s0MyM/KFeK5KxE2r2oBG8oxe83hVxQJn4DTi8PCBt0qwxvLO8jHEXDE24nanEp/E8ZjP738Kk27F0t8wjIu0xtK9+wGFgnohcrqr/stmmmEjoLGFeXi7HjD0Nxp5G6df7uP/vC5CafZx4/AhOP/0UCgoKQsps3ryVf/37FSq1HUMmXU2/1m0TaUKMJHcFK31CJv07m/7p4D3cQ/GsUQ3TxuKl0c0L2tL4nfr9oFsdTYp3vasxOESPBhfHJZQEeRJf09DuYF/04fCf+fQLvnHfuYt/P/sYt9/2w5jb6TDsWLqHyzMN+EJV9wOIyAvAJCDTBFbgw9+2Q2dGn3YtqsqGTSt5565HKGqTy3nnnsbQocewatU65s5/A2nZh2NPuZ6CZinbSEga0XxcHc3UFc/gq9dgLGOiPvDhDoMHHCSeHX/AUz/++sN51jB7V4iNhvNkAdPH5sCBI0p5eQWtWrWMUN7x2LF0fwX4sW99ayJQqqp7RGQHUCwiLTCmhFMxDExTQtJD1YsIfYaMo8+QcdRUV/Lvha9T8eQCWncdyvCzbiU3N8FbJkx4bSVSIAT/lr796MdHJ+a1Kqt1JytNyixcArxWxMHS2WtCfJNZ0agVxkZwmT5jzmXu3Fe55pqLY67LKdixdMcwKj0b2AJU4rNm9x14/g+wCqgHVmOEsU8JcUsPtWEi0KywBcedODPeW8SNedoQMKqPPz6mesxTD3CFlR87gsDKTYs5LV7L9eJZo3iXRSE+xszRb+LRtLR1B+qKQ+39Ovfoz/JPnsfr9ZKTY89s0ev1smjRO5xwwkTatm0dc1uaAhuW7gr8KEzZu4C7UtpAH0nXsJKJXUd0wWVCfFz5piMnLjg+opCN5McKsHQf42JNCSUBAsq8yxivf6nA0wyG0AsWXP6NgbHF9r2bCBrRB33H/ifyzjsfMG1a9PiWhw+X8pvf/oXW/U7h7Q+foXNb+PYV59OtW/xmOi6NSKTjNSKif33N2g4oz1tFS+8+SvP6WF5PFnZtmYJthoIpoYTyOYcRgbKBU2K+d7T7p4q29dupyOkU9iByJlJ0aBl538pvMIdYPcbYIR+9ahmzVu3hq9fGRKnBPivPvrthTSuS8WskRD20r9/K1/mDLK+rKmsW/I57f3dzxHpKSj7lkcdeZszZN9DSt5lUUVbKpo/n0oLDXHrJdAYPHhhSrn37Qaiqa9Nsg/gFllbRyrOXw3l9U9S0QKyER3BIKPPRkWBKKIFpsD5/Pd3Gt49Za/PfM904UWB1OryMUZeNaTA0hUaziGQKq2QhWk/7+i/CCiyAde8+y9UXjmbQoFDvJR6Ph3//+0XWfV7B6NOutZw61tVW8+nHL+E9spXvXnshAwb0bbjmCiz7xC2wcrWa1p6v0iawrIgUSMJPwO7RVNjQbANdx7WLWWA11VSwbf0OKnKKHCewci/Lb3gfj9aTTgyBtY2v80O1Hz811ZXs+PAR7vzljxvSvF4vr7zyBu9+WELPEdPpOTD6NNfr8fDhf+7m3t/c0LDz6Aos+2T0GlY0rB6Cg7MaDRjNgTIB1upaBp7Vh7IBsU8Js9WFcDxoy7Z89doQU0rmaVVmrFw1BtOssAWHq/IoLS2jTZtWvP76O7zxzgq6DTuT4gvOsX2vnNxcRp7+Q/7vnr/y27t/likeRx1DAhpWDa3rd3M4v1+q2pYUzEInb/wkBA+VudGt8TPFBsuJGlZR3WccyB8SPWOGkKN1tK3fwaH8yMfCDu79ks/eeYi8glZ0HjyFfsedEPc9d3y2jG55W7niipmuhhUDCRmO2hmZmhqzgMnzfN2ELXFxOh279GLASd+lS48BCWtGvYdMZPWb61m3LmXehI9K4vaHlfmiyhq73cxKk3KnhEcr9ntz154DYxZW+UsXsPLsu0PSR069mr898UpMdWU7CR7NcZbYitVbgyugsonUzMjyly5g8uwzmcyZvBu0DpqTk8PoM3/Mv5+0e3bZJW4Ny8XlaCGVi0eTZzceI6ornsFpZweezGzVJpJXGJdgEhJYTljDCsRd13SxIrKle7ysPPtuBn5nQIM9WtezV/Hma/lRSrlEIoE1LOdNCQ2c2GaXlJOCsew2bmPLo58zb1URF405kJFGs04jAQ3LgdqKA5vskg5SM4j5tamLxhxg7KpHU3KPbKOpfb2kGecFoTDiWrukntRMCbsynet+MwawfxjbJTyxCawGI1NFVRH1IuqhMWx4YFhwy/cNg1m4PI2vw4cj18a0kPpC8/rT87SSXK2jBfstygSXs7puymN5X4t8YQm8Hu5xKfBWkKMevF7z2kciGkH0BzOCKbGtunO0jlb1uwPSAl9aBZgPrtu8Qmr3mrleMSVZ5/f3olytI0frKPCWmdayGn8HlBVC8li97/LNdeS1yGPvvJGg3sA6Xev2uIkqsIrqPgt4b16czKOa9vVbfe8k6NE1vTd11NAvN0wZCJ9Pgjt8mM5l+Ro0aCashKkv7AMmwVVa3D8wPTJWD7BBjtZQk9Oa+tSFerNoSbjWBGPhwBEo9B6iLsfKQ6e12AisJpZBRE3Xwwx4Gu5ejd96jtaTo3U085ZGHhyt3qv1wFz779acyoVQv9V3Pbi8SzxEPZqTxra4uGQr21W1b1M3wglEFFguLi4umYRrOOri4uIYXIHl4uLiGFyB5eLi4hhcgeXi4uIYXIHl4uLiGP4fL1ZBvneFaEQAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"pysteps.visualization.plot_precip_field(R, map='cartopy', geodata=metadata, drawlonlatlines=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## pysteps using xarray"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"bom = xr.open_dataset(filename)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<xarray.Dataset>\n",
"Dimensions: (x: 512, y: 512)\n",
"Coordinates:\n",
" * x (x) float32 -128.0 -127.5 -127.0 -126.5 ... 126.5 127.0 127.5\n",
" * y (y) float32 128.0 127.5 127.0 126.5 ... -126.5 -127.0 -127.5\n",
"Data variables:\n",
" valid_time datetime64[ns] ...\n",
" start_time datetime64[ns] ...\n",
" proj int8 ...\n",
" precipitation (y, x) float32 ...\n",
"Attributes:\n",
" Conventions: CF-1.6\n",
" title: Bias Corrected Radar Accumulation\n",
" institution: Commonwealth of Australia, Bureau of Meteorology (ABN 92 6...\n",
" licence: http://www.bom.gov.au/other/copyright.shtml\n",
" source: rainfields 3.0.28 ho-rainfields 2018-03-23\n",
" station_id: 2\n",
" station_name: Melb"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bom"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<xarray.DataArray 'x' (x: 512)>\n",
"array([-128. , -127.5, -127. , ..., 126.5, 127. , 127.5], dtype=float32)\n",
"Coordinates:\n",
" * x (x) float32 -128.0 -127.5 -127.0 -126.5 ... 126.0 126.5 127.0 127.5\n",
"Attributes:\n",
" units: km\n",
" standard_name: projection_x_coordinate\n",
" valid_min: -128.0\n",
" valid_max: 127.5"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bom.x"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Simple Plot"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<matplotlib.collections.QuadMesh at 0x7f71cfe0c898>"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAZcAAAELCAYAAAAVwss1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztvXucJGV1///+zEz37s4uCyyLiIBCCJovGiVKEKNJQE0EXkbUeAGNoCYhMRglxvjTxIhizDfhF413cU24JShBIxGVgEAAgwZlQeSqERF0BcHlssvsLEzPzPn+UU/1VFdXd9fMVHVX95z361WvrvvzVHd1nTqX5xyZGY7jOI5TJGOD7oDjOI4zerhwcRzHcQrHhYvjOI5TOC5cHMdxnMJx4eI4juMUjgsXx3Ecp3BcuDiO4ziF48LFcRzHKRwXLo7jOE7hTAy6A/2grlW2mrWD7objOBXnER7aamZ7LuccLzpyrT3w4Fyufa+/6bFLzeyo5bRXVVaEcFnNWp6tFwy6G47jVJzL7Qt3L/ccDzw4x7cvfWKufcf3/sHG5bZXVVaEcHEcx+kXBswzP+huDBwXLo7jOAViGA3LZxYbZVy4OI7jFIxrLi5cHMdxCsUw5ryUiQsXx3GcopnHhYsLF8dxnAIxYM6FiwsXx3GconHNxYWL4zhOoRjQcJ+Lp39xHMcpEsOYyzn1QtJqSd+W9F1Jt0p6X8Y+kvRRSXdIuknSM0u5sEXimovjOE6RGMwVp7g8BjzfzKYk1YBrJP2nmV2b2Odo4KAwPRv4VPgcKK65OI7jFEg0Qj/f1PNcEVNhsRamtOg6Fjg37HstsJukvZd9IcvEhYvjOE6hiLmcU66zSeOSbgTuBy4zs2+ldtkH+ElieUtYN1BcuDiO4xSIAfOWbwI2StqcmE5qO5/ZnJkdAuwLHCbpaaldsqTUwCMK3OfiOI5TIAbM5H9v32pmh+Y6r9nDkq4CjgJuSWzaAuyXWN4XuCdvB8rCNRfHcZyCmTflmnohaU9Ju4X5NcALge+ldrsIOCFEjR0ObDOze4u+psXimovjOE6BRCP08/lTcrA3cI6kcSJl4AIz+4qkPwYwszOAi4FjgDuAaeANRTW+HFy4OI7jFIgh5goyCpnZTcCvZKw/IzFvwMmFNFggLlwcx3EKJo/Ja9Rx4eI4jlMghpix8UF3Y+AM3KEv6UxJ90u6JbFug6TLJP0gfO6e2PaukObg+5JeNJheO47jZBMNohzLNY0yVbi6s4lC65K8E7jCzA4CrgjLSDoYOA54ajjmk8HR5TiOUxmKHEQ5rAxcuJjZ14EHU6uPBc4J8+cAL02sP9/MHjOzHxFFRxzWl446juPkwEzM2ViuaZSp6tXtFcdph8/HhfW50xxIOike9drgsVI76ziOk2Qe5ZpGmWFz6OdOc2Bmm4BNAOu1YeCpEBzHWRlE41yq+t7eP6oqXO6TtLeZ3Ruye94f1lcyzYHjOE6MIRpW1Udr/6iqeL0IODHMnwh8KbH+OEmrJB1AVL/g2wPon+M4TkfmTLmmUWbg4lXS54AjiLKDbgFOBf4OuEDS7wM/Bl4JYGa3SroAuA2YBU42s7mBdNxxHCeDIkfoDzMDFy5mdnyHTS/osP8HgA+U1yPHcZzlMT/ikWB5GLhwcRzHGSXcoR/hwsVxHKdAjNH3p+TBhYvjOE6BmOHRYrhwcRzHKZjRHyCZBxcujuM4BWIw8qld8uDCxXEcp2Cq6tCX9MyM1duAu81stsi2XLg4juMUiKEqFwv7JPBM4CaidFpPC/N7SPpjM/taUQ1VU7w6juMMKUbk0M8zDYC7gF8xs0PN7FlEJZRvAV4InF5kQ12vTtJNOc7xczPLHPDoOI6z8qh0rZZfMrNb4wUzu03Sr5jZnVKxfe4lOseBY7psF1G+L8dxHIdQibK6Dv3vS/oUcH5YfjXwv5JWAY0iG+olXP7IzO7utoOkPymwP47jOENPhTWX1wN/ApxCpBxcA7ydSLAcWWRDXYWLmV3T6wR59nEcx1kpmKnKmsvBZvZB4IPxCkm/Y2ZfBqaKbCjXNyDpxZK+I+lBSdslPSJpe5EdcRzHGRUqXOb4M5J+OV6QdBzw7jIayhuu8GHg5cDNZuZVHR3HcToQFQsbH3Q3OvEK4AuSXgs8DzgB+O0yGsorOn8C3OKCxXEcpzuRQ1+5pl5I2k/SlZJul3SrpLdm7HOEpG2SbgzTezr2zexO4Djg34kEzW+b2balX21n8mou7wAulnQ18Fi80sw+VEanHMdxhpkCR+jPAn9uZjdI2gW4XtJlZnZbar//NrMXdzqJpJuJ5F7MBqJo4G9JwsyeXlSHY/IKlw8QOXtWA/WiO+E4jjMqFDlC38zuBe4N849Iuh3Yh6ga72LoKHjKIq9w2WBmpdjlHMdxRo35EpKfSNqfaET9tzI2P0fSd4F7gLcnB0oC9BpSUgZ5v4HLJblwcRzH6YEZzJlyTcBGSZsT00lZ55S0jshPcoqZpSN1bwCeZGbPAD4G/EfG8Tf06neefRZDXs3lZOAdkh4jGmwjwMxsfZGdcRzHGXYMMTufO1psq5kd2m0HSTUiwXKemX2xrb2EsDGziyV9UtJGM9ua2O3/9EjnJWDXvJ3OQy7hYma7tPWk6EQ0juM4I0JRI/TDc/afgds7BVBJejxwn5mZpMOILFIPpHb7pRzNzS2rsylyCRdJp5nZexLLY8C/AK8tsjOO4zjDThyKXBDPBV4H3CzpxrDuL4EnApjZGUQhxW+SNAvsBI5LDxsZhM8lr1nsiZLeZWb/NyQ4+zyRnc9xHMdpobj0LyG9VldJZWYfBz5eSIMFkvcbeAPwy5LeBXwZuNLM3ltarxzHcYaYeZRrGmV61XNJlsT8CPBp4BvA1ZKeaWauvTiO4yQwg0Z+h/7I0sss9sHU8kPAwWG9Ac8vo1OO4zjDSpXLHEt6OfD3wOOIzG2lRf72SrlfaH5/x3GclUCFTV6nA79jZreX3VBXn4uknikD8uzjOI6zUigycWUJ3NcPwQK9zWL/v6Sf0j1a4W+BrxTXJcdxnOGmwsXCNkv6N6JR/MkkxG2DM5dLL+FyH9Ar8/EPCuqL4zjO8DM4rSQP64FpWmu4GNBf4WJmRxTdoOM4zihjwGxFNRcze0O/2qrmN+A4jjOkVNnnImlfSRdKul/SfZL+XdK+ZbTlwsVxHKdgqipcgLOAi4AnENWF+XJYVziVFi6S7pJ0cyjduTms2yDpMkk/CJ+7D7qfjuM4MfE4l4oKlz3N7Cwzmw3T2cCeZTSUS7hImpT015I+E5YP6mMI8pFmdkgiLfU7gSvM7CDgirDsOI5TGSqc/mWrpN+TNB6m36M9g3Ih5NVcziIKW3tOWN4C/E0ZHcrBscA5Yf4c4KUD6ofjOE47Vmmz2BuBVwE/Iyqf/IqwrnDyZkU+0MxeLel4ADPb2ad6LgZ8TZIBnzazTcBeoa40ZnavpMf1oR+O4zi5MGB2vpoeBzP7MfCSfrSVV7jMSFpD9L0h6UASA3BK5Llmdk8QIJdJ+l7eA0O50JMAVjNZVv8cZ0UwVq+3LM/PzAyoJ9WnirnFJL3DzE6X9DHCczyJmb2l6DbzCpf3ApcA+0k6j6iATenx0mZ2T/i8X9KFwGHAfZL2DlrL3sD9HY7dBGwCWK8NbV+m4zj5SAsWpzdWMeECxClfNverwbxljr8m6XrgcKJUMG9N1WcuHElrgTEzeyTM/zZwGlEY3YnA34XPL5XZD8dZ6czPzLQIGNdaelO1xJVm9uUwO21mn09uk/TKMtrMW+b4CjN7AfDVjHVlsRdwYXDtTACfNbNLJF0HXCDp94EfA6V8MY7jLOACJT9mhZY5Lpp3EVUS7rVu2fQqFrYamAQ2hvEk8Te2nmgQTmmY2Z3AMzLWPwCUKdQcxymZsXp9hAWWmKuYQ1/S0cAxwD6SPprYtB6YLaPNXprLHwGnEAmS61kQLtuBT5TRIcdxnGGngj6Xe4j8LS8hepbHPAL8WRkN9kpc+RHgI5L+1Mw+VkYHHMcZLN38KR4ltnji3GJVwsy+C3xX0mfNrNGPNvM69D8m6WlEJY5XJ9afW1bHHMcpn6xIsG7RYfG25QqZkRZSFvldKsr+kv4v7c/yXyi6obwO/VOBI0KHLgaOBq4BXLg4zhCT9ZBPR4f1k1HxxVQtWizBWcCpwD8CRxINKSmls3m9Tq8gcqL/LNQDeAawqowOOY4zWMbqdbRuLarXUYaQKerhP1avZ5rdhn1cjRH5XPJMA2CNmV0ByMzuNrP3As8vo6G8gyh3mtm8pFlJ64kGLhauRjmOUwGe+otM77uOyStvw1KCpAjB0kl4jIrWEkWLVVZzeVTSGPADSW8GfgqUkkIrr+ayWdJuwGeIIg1uAL5dRoccJxe/dggTe3lauTKY/85tTP5oW+u6mZnSBEvajzMKAqYozUXSfpKulHS7pFslvTVjH0n6qKQ7JN0k6ZldTnkK0fCStwDPAn4POGGJl9mVvA79PwmzZ0i6BFhvZjeV0SHHycP4TXcwOzU16G6MLHO3fL85vxwz1VgwrSU1oNjUltaKRgWzQkORZ4E/N7MbJO0CXC/pMjO7LbHP0cBBYXo28KnwmcX+ZnYdMEVI4RVG6H+rqA7H5B7pI2kfSb8GPBHYTdJvFN0Zx8nLnAuWvrEcTSI+Nst3M8oUlXLfzO41sxvC/CNEOcL2Se12LHCuRVxL9Hzeu8Mp35Vz3bLJGy3298CrgduAubDagK+X0SnHcUYL7bYrAubv/3lTY0lrNKNEGaHIkvYHfoV2LWMf4CeJ5S1h3b2JYys3Qj/mpcBTzKwfafYdZ9GM1ev8513X8eRz38QB7/zmoLuz4ujkjB+r16Fewx6OfDijLFBiDDGfP/3LxriEe2BTyOjegqR1wL8Dp5jZ9vTmzG60Uq0R+gnuBGr0p4aL4ywa1ev89f1P46B/vKOc1zAnk17+mPmZGXhwIbw4aR6zmZmRFTaLUFy2Jkq4ZyKpRiRYzjOzL2bssgXYL7G8L5EwWejPwgj988ysL3+RvMJlGrhR0hUkBEwZBWYcZynMTU3x7UPG6VDexymJvP6Y5viVeg2x4MyPP0cnDJkwQr8Yh36o+PvPwO1m9qEOu10EvFnS+USO/G1xtd7EeS4ws1cB3wmVfVu7bPb0QjqcIK9wuShMjrOimdjrcdiOaQ8o6ML4unUt2kgsWFSvw0yjlLEzlaM4n8tzgdcBN0u6Maz7S6LAKszsDKKsKccAdxApAlmFHOMQ5hcX1rMe5A1FPqfsjjjOMGA7pmk868k89JTVPO4LtzP34EOD7lKlGKvXuwreUTSBZVGU5mJm19AjPYuZGXByj33uDZ93S3o8UVVfA64zs58V0tkUveq5XGBmr5J0M9l1lwtXpRynysxNTTF29Q3scfVC2KSzQCctZH5mJnPcw0hqLVQ3caWkPwDeA/wXkdD6mKTTzOzMotvqpbn0XZVyHGc0SeYNG1WhAmEQZcWKhSX4C+BXQtFFJO0BfBPor3BJqlJFN+w4zspjlIVKkqpqLkSRZY8klh+hdYxMYfQyiz1CF9eUma0vvEeO4zjDTnWFy0+Bb0n6ElEvjwW+LeltAF0i0hZNL81lFwBJpwE/A/6FyE73WmCXojrhOI4zOgwsnX4efhimmC+Fz8Kf53lDkV9kZslEaJ+S9C3g9KI75DiOM/RUVHMxs/f1q628wmVO0muB84m+tuPxYBnHcZx2is2KXAiSPmxmp0j6MtmRvy8pus28wuU1wEfCZMA3wjrHcZxSGcoIs4oJFyKXBsA/9KvBnsJF0jjwMjM7tg/9cRynJIYxxcrQljyumFnMzOJklZsJlYWh+XwvpWR9z2BsM5sjiihwHGfIGKvXGd+wO+Mbdkfr1jK+YffM2vVVpagKmH3Hck795wqiSpQxa4DLy2gor1nsG5I+DvwbsCNeGRexcRynejQFyEwD6jUANDnZXDfGgqmp6qanqvYrE6OKZrGY1WbWzM9jZlOSJrsdsFTyCpdfC5+nJdYZ8Pxiu+M4TlE0H8gzM4yvWweAMY1qNajV0NpJ1GhgUzuaqe+rLmSGhQoPotwh6ZmxYiDpWcDOMhrKm7jyyDIadxynPyxUf4w0GGs0IiEDaN3aSLtximO+sprLKcDnJcX1XvYmqjJcOHnLHO8KnAr8Rlh1NXCamW0ro1OO4xRL0/w1BcTVIRvtAsU1lmJor5hSDczsOkm/BDyFaED898ws881CUp4yKw+a2euzNuQ1i50J3AK8Kiy/DjgLeHnO4x1n6Eg6vUfpoWszM5CoY99c5xTD4Jz1PQn+lbcBTzKzP5R0kKSnmNlXMnb/P8AfdDsd8IlOG/MKlwPN7HcTy+9LFK5xHGcISQuWURKgg0VVduifBVwPPCcsbwE+D2QJl78ys6u7nUxSxxH/efNC75T0vMQJn0tJTiDHqQpxGOwoPHTj8GPFFSGJhMrc1NTIXGOlqG4o8oFmdjrQADCznXQoRmZmF/Q6Wbd98moubwLOCb4XAQ8CJ+Y81nGcATKWEijgWkrpVNQsBsxIWkPooaQDgceydgwDLP8A2Be4xMy+kdj2bjP7m24N5Y0WuxF4hqT1YXl7nuMcx6kG3UoPOwVjVDla7FTgEmA/SecBzwVe32HfTxMNuPw28FFJV5vZ28K2lwNdhUsus5ikXSV9iKg05n9J+mDQYgaCpKMkfV/SHZLeOah+OM4w4FpK/5Hlm/raJ0nA94gEw+uBzwGHmtlVHQ45zMxeY2YfBp4NrJP0RUmr6GBKS5LX53ImUcWyV4VpO5FjqO8EVe0TwNHAwcDxkg4eRF+clcEwpUtxKkIFfS5mZsB/mNkDZvZVM/uKmW3tckg9ceysmZ0E3EikZKzr1V5e4XKgmZ1qZneG6X3AL+Q8tmgOA+4I/ZghKgPguc+c0sjr8HYB5AwB10r61Zz7bpZ0VHKFmZ1GpFjs3+vgYYwW24fWms9bwroWJJ0kabOkzY1sf5XjFEqVzE8u6AZLFc1igSOJBMwPJd0k6WZJN2XtaGa/Z2aXZKz/JzOr9WpoKdFiAA/R2QlUNlm2vqziN5uATQDrtaG6sRuOUyBj9Tpat5a5Bx8adFdWNtUd53L0cg6WtCmYx3oyjNFiW4D9Esv7Avd02NdxVgxxyLELlgFjwPygO5GNmd0t6ZnA8wiFHxeZ3T6vSS13tNjfStrNzLab2XZJu0vqGoZWItcBB0k6QFIdOA7IkwPHKYm4RogzWObDoEhn8BRlFpN0pqT7Jd3SYfsRkrZJujFM7+lxvvcA5wB7ABuBsyS9exGX9rO8O+b1uRxtZg/HC2b2EHDMIjpUGGY2C7wZuBS4HbjAzG4dRF+cCNVqjO2+26C74TjVobhosbOBo3rs899mdkiYTuux7/HAr4YArVOBw4HXdjtA0qGSLpR0A7BvNz9Nkrw+l3FJq8zssdDYGkoqjZkHM7sYuHhQ7TutWKOxUITKcZzCwozN7OuS9i/mbADcBawGHg3Lq4Af9jjmPOAvgJtZhMEvr3D5V+AKSWcRfW1vJFKtHAdmGsz94gbY8tNB98RxBs4AIsGeI+m7RL7nt/ew5DwG3CrpMqJn+W8B10j6KICZvSXjmJ+b2aJdD3kd+qcHNeiFRNFa7zezSxfbmDOa2MwMjfWrcr+pOM7Ikz/9y0ZJmxPLm0Kka15uIEqfPyXpGOA/gIO67H9hmGKuytHGqZL+CbiCRB4yM/tit4NyPw9CvHNbzDOApP8xs+dkbXNGn/mZGVbfuZXZQXfEcSrCIjSXrWZ26FLbSUbumtnFkj4paWOnkfdmthSL0xuAXwJqLJjFDChGuPRgdUHncYaU2TvvGnQXHKc69MksJunxwH1mZpIOIwrSeqDgZp5hZr+82IOKEi4+SNFxHAegQJ+LpM8BRxCZz7YQZTWuAZjZGcArgDdJmiXKmnJcyCFWJNdKOtjMblvMQW4mdxzHKZriosWO77H948DHi2mtI88DTpT0IyKfi6Km7endDipKuFQ214HjOMUQD5Ttdw61TgN0q5TLrY2K2XIkfZkuvTKzl3Q5vNc4m0xyCRdJbwbOC4Mns3jdUhofFGP1erVvTMepGIPKwBC3m66kGW+r6v94QEkpu/EP4fPlwOOJhpdANKjyrm4HmtndS2kwr+byeOC6MELzTODSpF3PzDJTE1QRT1PiOIujSv8ZJfpiFRUsQOU0FzO7GkDS+83sNxKbvizp62W0mSv9i5m9myh2+p+JsiH/IOQbO7CMTpVF/KZT1bcdx6kS6SJp/frvxO3mEWqVLOSWM6/YgLSbPSU1a3FJOgDYs4yGFjPOxST9jChx2SywO/AFSZeZ2TvK6FyRVFmFdpwqUhXfinIKmUr9vyumuST4M+AqSXeG5f2BPyqjobw+l7cAJwJbgX8C/sLMGpLGgB8AlRculbrxHMdpoZv2kTR/ZfleKklFhYuZXSLpIKJBkQDfi3NGFk1ezWUj8PK0Y8fM5iW9WNLuXZz9juM4hVB5oUIUOltBhz4AkiaBtxGljPlDSQdJeoqZfaXotvL6XN7TKWLAzG4nyjnjDICkfbpytmfHWYkYaD7fNADOAmaAOF3XFqCU2lw+zqVipB2ovfZNR88MaiyC4yyHpd6veV+osvYr9T9SUc0FONDMXi3peAAz2ymplOe3p3/pI+kbvOg/lOr1vpoNXJA5gya+9zr9JwamzVf3iTgT6nEZQIj4HajPxSmYXg/k9Pb0YDIA6jWYaQzMDu1CxakSvaLNkv+Tsu/dqvpcgPcSZbffT9J5wHOJsh4XjpvF+sRizF3zCfMWhD9HvYZqNajVIqFSm4DGLMw0YHq6RcjkMafl2c9xRoF+CpWFRvvTzGIxs69Jup6ovLGAt3ZKz79ccgsXSePAXsljzOzHYfYFBfdrxdJRsExOYrutY34yvIVNjAMwtnOGsW017P6t+X00QeMZwwWMM7rMD8oHaQNz1vdE0hVm9gLgqxnrCiXvOJc/JUr1fB+txWKeDmBmDxbdsVEira7nGfAVq/JaOwm1GrbbOhp7rGV2zThWE/MTYmzWWPUAjG2bhnqt63ljwaK1Ua17A3DB4owIne77gb08VUxzkbQamCRK3b87C9am9cATymgzr+byVuApZlZ0EZqRJyuiK17f1a8Sm8F2W8/crmtorF/FY7tP8OjuUfT4/ASsftgY3znR80ccq9fRurXR+QBrNGCmgep1114cpwQq6HP5I+AUIkFyPQvCZTvwiTIazCtcfgJsK6MDo0JaO+mUsiIZ0ZU2f9nMTKtgWTvJ/GS9KVhm14iJR43Z1WJsFiZ2zjPWmIfaRLR/hrBICjdrNKKVM40Cr9xxnDYqJlzM7CPARyT9qZl9rB9t5hUudxLlo/kqibA1M/tQKb0aMlq0k3qt5/6doh/aHPe1CebX1JmfiI6oPzLP7JoxJh41JnYu3L1zu65hvDHbPG9awNjMTNMEltaiXGtxnIIxKidcYszsY5KeBhxMojy9mZ1bdFt5hcuPw1QP01Cz2CR33WLl00KlKRjidTMNaKQ0haR5qgtWi5z29e0NVj34KACN9auYnxBza8YYmzUa6yeYr41hE+NMPDARlYgLjsysLLbJlAwuWByneER1w2clnUpUNvlg4GLgaOAaYDDCxczeFzq2S7RoU0V3pJ8sR7C0mbuSQgUWBMvaNWF5AhpByMRaTZbASWCNBqrVUGOOsZ3BhDY9g9XGqQFza2qMzRpza8aYnxCaUCRgauPRjR1MbFlC1AWK45RPVaPFgFcAzwC+Y2ZvkLQXUTLiwskbLfY04F+ADWF5K3CCmd1aRqeqSm7T19o1WG08etg35qA2jmrtX7XI0F5if8gk0JhlbDqsb8yiHTsZbwShRY25NWOM74zu4qZwCeNgkhoMuFBxqkPl0uOXQUXNYsDOkHB4VtJ64H7gF3odtBTymsU2AW8zsysBJB0BfAb4tTI6VSXiWPleNSVibaN5XDweJQiY+DMSNhOtTvUw3zLSfnoaNRImtkYjmt+xE03WmWjMMb5znMb6VUxMzwKRdpMk2WePCnMGzYpKrFpd4bJZ0m5Ez+/rgSng22U0lFe4rI0FC4CZXSVpbRkdqirNSC5oj7bqoMnEAx2BSKgQhM10FBMRhwSnRxCP1evReoK5LRYs4TMWIvOTdeoPTGO18TbB0mauwwWMMzhWmmCpYCgyAGb2J2H2DEmXAOvN7KYy2sodLSbpr4lMYwC/B/yojA5VkfiBbxkRV0CrPyUlaDQ719RYYgHT9LkEwZLlFxkDxtZOLggWWBAwO6J2xqZpFSyNWYg1qDBYspd/Jw/p1DUrwqzhtLGYFEbdWBH3TsWEi6RndttmZjcU3WZe4fJG4H3AF4lcBV+npGRnVaVl7EhqW3IciYJWk9QkYk0lftBbF8GSxBoNNBmERHxcmI9NZkoGCUAkVGoTC/6eHdNYo3dbncjKLqB6nfGEsF0RD4sCGFb/V/ybLyZB6rBea1FU0KH/wS7bDHh+0Q3mjRZ7CHhL0Y0PK5lpvmO/CdOIyYVQxHj99HTrcpeHfTNVSy1lbosFSD3hh0kSJ7QkmOEe3s78Qw+39HkxZJkykrnJlnrelUjyuxxGzS9v8sesl5Ekw3bdS6VqZjEzO7LfbXYVLpI+bGanSPoyGYqemb2ktJ4NAWlzGRANVkyaySDTr9KJluSScQbkHdPYjkg42cxMFAlGwh8DCcf/7IK2MrUjs6085o22FP8ZfqWiHxQr5W132K9vMcJxaGreF0mBgyglnQm8GLjfzJ6WsV3AR4BjgGng9d1MXJJOyFo/iEGUsY/lH4pueBRJCptYc1lsGnygmQW5Wa+lsaDtNAVPfP6wTdBicpubyh6KtBjHakvyzERb3c6/VEZdsCR9VcPyNp9ZQwg6jqFK0iu6cuQpTnM5G/g4nQc5Hg0cFKZnA58Kn5341cT8aqKM9jd0Of+S6SpczOz6MHtIyE3TRNJbgauL7lA493uBPwR+Hlb9pZldHLa9C/h9YA54i5ldWkYfFst8Woh0SVCZJGuQZtMclnLGt6SHSZHUbIp8WGntZOTHSfSljLfQqj5giyR5jVV/o89MdgrNhKe9+m0Z17pSEMWZxczs65L277Jw2LK2AAAgAElEQVTLscC5ZmbAtZJ2k7S3md3b4Xx/2tJXaVcWlIhCyevQP5FI9Ury+ox1RfKPZtaiMUk6GDgOeCpRds/LJT3ZzOZK7EdH0nZ0aH9ILvmhGT/IgyYCC+n3W/bJ4cNZcn/qtWaAQLOdDqY2Jz/ziVxvVWR8w+7tQiUmaNOdMmpn+SOTgmYY/U1Lon8+l32IEgvHbAnrMoVLBtNEWk/h9PK5HA+8BjhA0kWJTbsAg0i/fyxwvpk9BvxI0h3AYcD/5D1BUeaXTqn0iyAr51hTsCTylcX7FWmiappC1q2NTHO1CXh4O7Zj2hNdrgCaY6w6jeVKre8kLLLMfysGA83nli4bJW1OLG8ys02LaC0rjVnHxlP+8zGiHGMXLKK93PTSXL5JJAE30hrK9ghQysCbBG8OzqfNwJ+HiLV9gGsT+8RSOhdF3exZI/a75fPqRvpNL+mvaZIeQ5PDt7IUMm3sLlhWHOkXJQWBk6bXPZ8ZVcnK0F4WYRbbamaHLqOpLcB+ieV9gXu67J+0Bs0Cd5vZlmW035FePpe7gbslvRa4x8weBZC0hugi7lpqw5IuBx6fsemviJxS7yeSsO8nEmxvZBFSWtJJwEkAq5lsri/spu7wJrfU83cfR1NbGENTolM9Jk7R7yawlUdWAbuOg4dZ+ktV3gGZWS+EQ3FP9s8sdhHRi/j5RI78bZ38LQBmdjVAyCs2EeY3lFFNOK/P5QJa84jNAZ+nNfJgUZjZC/PsJ+kzwFfCYm4pHVTLTQDrtcGgmJuyTWuJI7oKvOHTxb5sagdat7Y5VsamdhTWVqd2HQdazVstKZASL1edBEy6GF6apfxnhuUeLcqhL+lzRCnyN0raQlRuvgZgZmcQpc0/BriDyH/SdXB7eOl+P7CTqGR9PLKh8OSVeYXLhJk1f1Uzm5FUmkE1Fe3wMuCWMH8R8FlJHyJy6B9ESUnXssg0q5Vc1bEZ3jy1o/knHpY/mDNaZNUuMmgTMN2OaZL43+TVeobqvi8uWuz4HtsNOHkRp/wL4KlmtnVZHctBXuHyc0kvMbOLACQdC5TZudMlHUL0E91FVP8ZM7tV0gXAbUT2wpP7FSmWrm+fRVk3fzq82SmP5O+8kgV5ptaRURAvfu2N59tI/FdUq3UskJclYIrKZdZ3Kpy4EvghkYZTOnmFyx8D50n6BNG9tAXIHOlZBGb2ui7bPgB8oKy2s8gjWJzhpyU7QpwjbsQHd2bRzZzVUmm1F6n/SnKsTK+MFUMrWAjjXKqXWyzmXcA3JX2L1pL1haf3yptb7IfA4ZLWATKzR4ruSFVp8bFkvYUlHkTOcJL58lAPSUEL9qcNE2nfopK1hbJIB7kk/xeJ/0neVEhDjVVWdfk08F/AzUQ+l9LIW4lyL+BvgSeY2dFhMONzzOyfy+zcoOkkWGChOFgnNd+pPpnBGUlyZK4eVbK+l8yqqR2iJiHlsO8SqtyNYf3uK2wWmzWzt/WjoXTUayfOBi4lcqID/C9wShkdqiwzjbY/UHL0vDM8jNXr2dVFU9VBFyNY4nOO5MDB+N5PT8ltKXppe72+1/lhFuq2iKn/XCnpJEl7S9oQT2U0lNfnstHMLgh5vTCzWUkDSbnSL5Lhl0kynZYFtZfM2bSUMQND+2fsM51Gj8fje5L7dSPt/B8V4mJ13XKC5ckdtpLvxwr7XF4TPt+VWDfQUOQdkvYInUDS4cC2ojtTJTpFrmSNXl7OgyXtRI7/qL1KEqdH06fDQFfyHzsPye98KaavTrmzRuF7X3Tm7JTfcaWaEpNUVbiY2QH9aiuvcHkb0RiTAyV9A9gTeEVpvaognQaHFSZYoM1+3S08Mx25lqxIuZQR0yuFdIG3pT4IR/27zfRFJQRImy9yhQc/tGBU1qFfpXouccM3SPpN4ClElqHvm9mKdTakBc1iTVi9nMjp/WOSf+iWkNDkm2NtIYIt7VAb9QfiYhimEs0De1FI35vdwvBTgmWlv9xU2KFfjXoukp5vZv8l6eWpTU+WZMCDwDWDSnlfBRbzBxpft66tQmXLZ6CjPTspVOLzhLLG1CaiKpQsDGxTKvxzpf/h01T9uxhEcECbNg3Z9YMarcEPaar+3ZZORYVLleq5/CZRTPTvdNi+B/Bu4LeK7NQoMtZpAGYHU0KWdtNJsFhtHDXmFgQNC4EHydHTbrJopcrCth/+s+T1dzS3JgVLYl5kl4bIOvdKo8hiYX1gMPVczOzU8NkxGZqkkR7rUgSdwl7TD/vkH70lSSCpkdFJjade65yiv5GtGa1EqpZdt1ep4zL6lvUdZAmVNoFSry28uATtGPIJGFiBWoxZlX0ulannEndoV6JsnL8RVl0NnGZm28zs98vo2KiRrtPSLcKomawyJWCS1SlbHgDxYLb4AZAsJpazuNNyqPpDJB3mDYPra1aUXz/60snE1lGoQLMKqdXGF/aHFgHTdq4ElnhZquq9URZVjRajKvVcEpxJlJn4VWH5dcBZQNoX42SQFhaLfshlVAZsLsd+mEajqdnEqfmz0m0UTZUHDWZF9Q1aWymremm3NoHMdmPBkixp3NSQ167BauNNwaJG5FaNTLDZwiVt9m2aZleYYIFKm8V+DNybrM0laX8zu6vohvIKlwPN7HcTy++TdGPRnRllkgImXs5DrgR/A3qIVlmwQHslxEELFsh+0JalTXbSJNJ+lVioaHIS6jVschXzkwmTbGOu6deLhUwLOZK5rijtxYD8ZY77zecpuDZXJ/IKl52Snmdm1wBIei5RsRmnJPL+EecTo8oX8wcu+s/erwfHUvo9P2DzTKbPLVCGWayrCQzaNYxEoIhNrloIEEkf35iLTGLJiqg5SA4MhuqaTwulsrKlf7W5FpNy/9zgewF4CDixjA6tBEqv+9KnPvT7IbEcTWkQgiUd1tuS6HRQQRZZgiVen4w2DNpKzNj0DOyI3iebZte8bYVr7VRQbBSFTYXNYn2rzdVTuEgaA55iZs8IdZcxs+1ldMbpH3lNRd0e6P1+KCynvX72NSlUmg/v2C8GbSPde6X6WSxtASEdzFZttVkas6gxi4ITv2kGa8zCjumlZQBPpd1Xui8Zg33T1zKUVDRajIXaXB8Py6XV5uopXMxsXtKbgQtcqIwOi/nTdjLprCgzRw7anOczjWgw6+Rka/Re2wM2O8VPYd9rIvijjbTAiLWs2CwWzGBtkYe9/CwZ+2YOxoSOUZQwpL4aq260WD9rc+VNuX+ZpLdL2q/sNM2jzLClEW/TWuq1XM7bskimta9SMEHcn/j3nZuaWnCeJyKxgMwHrVLXVMR9Eh/fWlMlo2REmplGJFDiaZH+lZYU/OF+aQYLrJ2MNKVkmeSUwFHQ+pJT1X7vXkSDKC3X1Pe+SX8raTczmzKzRyTtLulvymgrr8/ljeHz5MS6UtI0OxUlKVRSWXDLfrusuo2+Y3/qwceSLiwXNJrkfn3JoNDJ3wJNzUpJ8103evmMslIVAaxd0zpOptFAayczBV38ncVpjIo2H5ZKRTUX4Ggz+8t4wcweknQMUaaVQsmbuLJvaZqd6tBS96RL1cGl/uHzCqWyHyhlpVqxqR3R23f88EwVI0tThv+lzeeS+B2zHujxuqTgifuep4YL0J5YNaY20bYf0CbMksf2ygJQVQahleRkXNIqM3sMonEuwKoyGso7Qn818CfA84g0lv8GzogH4jijTdMxXGAN9CqYOTqlQykqkm583ToAbMd05/xxWYlMC6At/LnXubv5U3pUmuxUErklLVHynMlUMvUFIRJtS0WvNWpoptHM9j0+U69+vZjBVZnMw78CV0g6i6iXb6SEjMiQ3yx2LvAI8LGwfDxRJs1XltEppzqkB3+mtw0r3YRbOqHjUq9zbmqqd/uxJlGQmTGdlSD5UgB09Zml0w1laSxpupWKNrKd+E3iTN5pX15Sw2kKomHSXgxVdBClmZ0u6SbghUQy/f1mdmkZbeUVLk8xs2cklq+U9N0yOuT0n17RSemR7unjuh3biSqMls+7T2lmsy4a4VLKXqevq5vpqukDykoRlJFuKPMcPUxjLaRylAELmbwbsy0CpSWXWRy1Vov9Uo3heKmprlkMM7sEuASiAfGSPmFmJ/c4bNHkFS7fkXS4mV0bOvRs4BtFd2aYGMoQyQTdMuTmva7kw6UKYcl5TW2d8nt1EqJl0GyrS7tLcfC3mNqC1pEemd8SXJBRorjruRP7tSVWTRI0sqZ5rFGLcpIlhExTiKRymHXKEBD3ufL/vYJDkSUdBXwEGAf+ycz+LrX9COBLwI/Cqi+a2WldzncIkfXp1eGYLxbX2wXyCpdnAydI+nFYfiJwu6SbATOzp5fRuSpTZsqOFkf6MtrK85BMZufNe1zWAyU5+rpTf+Nt6X2SQm2pYzyyEjS2ED9EO4waT9KppHUZD7T5lHDLSm6ZV2tJ17PPJZxylibO+p7ynLvl+06MoWkTMIH0cprK+1tiCtJcJI0DnyCqmbUFuE7SRWZ2W2rX/zazF3c5z5OB44iEygPAvxGNczmykI5mkFe4HFVWB5zu5qasZej+oMssjRzo5IjtauKI6VHmNj7HeJcHUdy3bppTcn4x5iCtW9t9xy7h1HlMgmUnwGwRLKn+5aZTVFdaa1lkvfv0Q73TC0jPiLKgvaTTyzSPSWgszfmZBjY9jU3tGA7BAkU69A8D7jCzOwEknQ8cC6SFSy++RxSI9Ttmdkc4158V1ssM8oYi311mJ1Yyee3kWSOXs2iL3Ek9oDLPnw4d7SREko7WHTujsRE7EnmmwnmaYxPSpFOBpAYYtphrwv69QnPj0tGq1WDtJDYZRVVq+rFmypK0UznL17DYpJ/06NdSWI5gmZ+ZaZom237jVAbkxQiWbkI3q6RBR8J4n2ZZiJAgsxPJJJlDJ1goNBR5H+AnieUtRJakNM8JfvB7gLeb2a2p7b9LpLlcKekS4HxorzNYJHk1FycHS3nopP+knchjluiU9bYjqYFu8YM5XcMjvS7642cLFmoL9WUykzRmjIFoPvxTgiBODZIlYOK3fK2dbAqV+ck6NjGOZucYb8w1ncBJ0mMmlmTuKYGWeyA8/Jf6MG1JrZ+mAMGS3N4pxUdb1BmprASJ+jBtx6YEy9BVUjVgLrdw2Shpc2J5k5ltSixnCYD0yW8AnmRmU2FA5H+QKl1sZhcCF0paC7wU+DNgL0mfAi40s6/l7XBeXLiUQFkOx24PwkUJlmQKklBtcH6y3hal0/HNsjHbFCwtTt1QZCpKHbLI7L/x/uk67QApP02LYNltPbN7rGO+NsZ8bYyxxjwTO2eafeioFS2B9Pdf1O+crPVSRHLO5Pni3yZ5/UUK1LwvR4smTjuzTEE7CMSiUrtsNbNDu2zfAuyXWN6XSDtpksz5aGYXS/qkpI1m1pbt2Mx2AOcRJa/cQDSc5J2AC5cqk/yj5X3wLPWP2c300ZUMbSVdcXBseia7lG1sosh4+1Wc5iPOoBsGDqYfeHE4aWYqkrWTXXNZJX0rmpzEdlvH3Po1rYJl+6OMbZtuy+KbpUUt9gGb1giS6eOrlLG5JTCgy2j8ovqUxxfTJNEP0R52rMbcgsl1EX2uXIqg4sxi1wEHSToA+CmRaes1yR0kPR64z8xM0mFEyv4DvbtoDwKfDlPhuHApkaI1mExzRw8ne0dS2W/VmGvxU2SdKx222/KHjjWODMES7x+dJxyXdvzWahBrQzta64WkrzeumAgwtnMGzY4zvjOqOdK8hh3ZNUe6CZRuD6jmA3vd2gW/UvLYDufM+v3LjkBLm9jaIuV6UGifUulm2hJVph34ScGSQ2tZanRd6RQkXMxsNmSlv5QoFPlMM7tV0h+H7WcArwDeJGmWqIjjcWaDH2jjwqVg0g/dxQqYRb1N11v9FvGft9MYhhaSNUV27OyokcRkXUPyWm16GqbJ9UBIbxuv11uTOnaIdmrx1cw00EzrA96mp5mf2tHSRjeBkfWW3cv0GD8cNTkJkywkfMw5XiQrErAvD8P4XllqFNpySI/aT26Lx7SkNeUcfexW4XOgGIUmrjSzi4GLU+vOSMx/HPh4+rhB48KlBPKYCYq0UzcfeBnRVm3Eo6/jt8n7Wx/Gi6Xp1J3JP5I8Tac0KU0zWJyqPRA7ebulV0n2bzF0EoxxVBrQEpUGiTfuWqMlj1i3QZllDkBtecFJFiyLqdXaMxEX7NvomEttptHqV5t+rPU+zWEOyxzTlMhRVoXsyRVOXNk3XLiUSLe30k6O0LzjV7LIelh0pMAHSRlO3eTDsbVa4uLHaHQi/RLQbfBnk1qtOco8OZK8mb03kT6/l1Dp1Jfl0rESJrRoLm1RfQXQywzYJBWa3FwXE/qZFSXYQkp4ViM9jMF8dXPu94u8xcIKR9IrJd0qaV7Soalt75J0h6TvS3pRYv2zJN0ctn1UUqlx2sslTyhnespzvvSD1ZJVDmca2NQO5qamMs+ft628lDm4sPkgboRAgunphQf4urWMr1tXmFDrZcabm5paCDZozC44n2Pih1yGWU+xr2bd2mi/tQs+o4GEQXdIdQ+tv+dySZ8rvlYLfj0Lv2uLoKllf4/pc8XRgppc0Gptx3QubbZ0jMjnkmcaYQapudwCvJxUpIKkg4kiIp4KPAG4XNKTzWwO+BRwEnAtkQ3yKOA/+9npqmAzM5HtOmVD7+efaykPobyCqMXc1iXMup9vqTYzs6CZ1MaZXb+ascY8Y0S/RWYgQj3ho0n6yCD63YKJryhTTtt4GVKZieOgiYzyxvGbWpkCryVZZ5ZmFfoW38eZ2l+ciSEWKg9vq4ZQSeKKy+CEi5ndDpChfBwLnB+K2fxI0h3AYZLuAtab2f+E484lGgy0ooRLmwlqGaG1yyFv6Gc3n1MuARMc7HGUFix+lHnh1CaYX1NnrDHP+Pad6OcPtdVsSY4viccSWW2cufVrgCjKbWx6BjUaLYNFkyxV2KTHPFmj0RQqyfUtARKNRktAR54Xh9wmsNT+Y6TMg3EYe2JsS6dzJYMqaDSYf+jhCpjB2nGfSzV9LvsQaSYxW8K6RphPr89E0klEWg6rmSy+lwMk+SfNWt+vPnRKS59FVo6zXgIm7bjtmCK+TzQfbDt2MvHw9o5mmLHY3zE5CWvXMD9ZZ2aPSRq7jDO+c56JnXOMb2+PkEonEU2n3YfOD/Skj6XFRJoWeGkfRWoEfM/xUxmO//SYmvRv1CnAJR391/T/1GuMr1vX8Te2RqP66WBcuJQrXCRdDjw+Y9NfmdmXOh2Wsa4tgjGxPpOQQmETwHptGMlfejlZhMtgsSO/8/S5+bDL0ND6fc1zU1OwBPPL7PrVNHYZp/bIHLXtjzG+9ZHW8URZ+d/SJjSyBXJy7A20mrvafHPBlNot5LzZPrTnf0sVNEv2VQm/R7LtrN+orQInNP0/vTIpVM78lYUZzLldrFThYmYvXMJhndIdbAnz6fUrmqoIlV4sVtNIpzCpQv6vRVOvMbdrZAab/Ol0u1CB7HDxdPhwKtIrqbEmB3V2PSfhXulhvswSLMl+NbNfJ01rgaRg6RnlFYckJxz56cGTbX0fJlxzqaRZ7CLgs5I+ROTQPwj4tpnNSXpE0uHAt4ATWCi77PSZNpNVwdpEUisbirfVBFq3Fu22Kza5ivF7HmTs4W3AErKwx76I+LwzDYzpVhU+VTVyOQ/lFi0k2XZtIir0FdpLR5nFZqpke10Hr5J6IUlGO6auYeiESowLl8EJF0kvIxIOewJflXSjmb0opDa4gKhewSxwcogUA3gTcDawhsiRv6Kc+WWTd1xG8s9fZrqNYXuwJEeM2/1bo2qNie3dNLyOqX1qrX9RMRkN1kyZ1LI0uyULljBItDXvV6I6ZML5nhQssBCokaftXvnehu33b2LAvAuXQUaLXQhc2GHbB4APZKzfDDyt5K6tSLqm0sgwkVQmh1OF6GR2yspSHJPpF0n5IxYSi9ab42sU5+AK+4ilP5ibWmitBrutZ27XNdhESGQ6O5dZcjhKUtrqh1n0/dBBMA7/fWVg7nOpollsxVBUAsOsCoaLLo+bODads6wt9YqTm7Q2mCfyK5kiZX6yzvyaELQRMj9rdo6x2jiafiwylaV+kyXdS3EmhMZs5BuqtT8arDaOhdIM49t2Nq+nyJQ/I4ObxVy4DIpupYx7/VlbwjnjSKFkLZWZBmo0UM4Ryy0D22AhPcjaRAh3lygkJz95xvaMp+4Nzc41NYn52hjjs1GtHdUmot+6Vlt22hOb2oFN7Wh31idKM8TlDSa2Pwo7dmKNpbU5/JpJDwyPFsOFy8BDeTuN/8hzXDOhYxhLEVFvps8XMJ5Tk+mYCXjtZEv+qbTj1unOUr6nZiaAmQZj0zPNiLP52hiNXcaZn6gzNmtRBNqOnctOe9LSx5mZxPiaRqj0Gd1fmp1jYucMYz/fhj28zV8yuuGay8oWLstJjb9c0gMhmyO6M3Ixpdc3hUqHYl8ANrkKkRgc9PA2xqZa207PJ4lHxmf5CVyw9IGZBkxGZrFH91rDzC7jPLqbsAlYe+98JFjuvq+UEepJ39E4wI5pxn66cD/O+u/fg9HPG5aHFS1cqkAyxUlMelxHOoIoWahK6WiiUJ54bHomspvHZrJaDcJ5l1JjJtmnfgviFUujgRpzrHpghtr2MXa5az5KG3PPVuYfergvD/mR942UgeFZkVnhwmW+Ig/L9Ijr7AiihZHqkdkkaDqhQFgziWVtojW6J5lWY4lUsiDTCBNrDmMzM2jHNONEGkRs2vTH1hDgmsvKFi4xVXoL72imSm+LBWN6fERGLijIH6aaTmueNw2JUzzdRtQ7FceFiwuXYaFb0TFIJLHskC5lPodA6CZYgMGUyHWcYcMMm8sYG7TCcOEyImRpNV336bA9M51/QqgMKmmk4yyXZJbp0u9fH6HvwsXJJunnSab3ABcszvAxvmF37El7o+kZtGMn81t+Wm6DbhZz4eK0kgxyGIk8T86KZqxeZ+yAJzK36xoa61ex+p6tzN53f7NeTCn3tZlHi9Feb8pxAJrJB/MmIRw0RdZ+d0YLq43z8C+tY8cToowD2054TmbJ7GIbtXzTCOOaizPUuEBxOjG+YXcav3wAOx9XZ839s6x68FEeeP4T2fDFm0v2u7hDH1y4OI4zwtR/8iC1B1bR2GMt0/tMsselP2S27IGhnnIfcLOYM8R0K0jlOHMPPgS1Caw2Tv1n21lz4beYve/+/jRu8/mmHEg6StL3Jd0h6Z0Z2yXpo2H7TZKeWfj1LAHXXJyhZNAJR53qM7HX4wAYC078fmGAFaS5SBoHPgH8FlGp9+skXWRmtyV2O5qoYu9BwLOBT4XPgeLCxRlKXKg43RjfsHuULbqPQqWJFVos7DDgDjO7E0DS+cCxRJV6Y44FzjUzA66VtJukvc3s3qI6sRRcuDiOM3LMPfjQQNsvSnMB9gF+kljeQrtWkrXPPoALl7J5hIe2Xm5fuLvPzW4Etva5zTIZteuB0bsmv57l86TlnuARHrr08vkLNubcfbWkzYnlTWa2KbGs9AFAWnLl2afvrAjhYmZ79rtNSZvN7NB+t1sWo3Y9MHrX5NdTDczsqAJPtwXYL7G8L3DPEvbpOx4t5jiOU12uAw6SdICkOnAccFFqn4uAE0LU2OHAtkH7W2CFaC6O4zjDiJnNSnozcClRWZ8zzexWSX8ctp8BXAwcA9wBTANvGFR/k7hwKY9NvXcZKkbtemD0rsmvZwQxs4uJBEhy3RmJeQNO7ne/eiEb8fw2juM4Tv9xn4vjOI5TOC5cCkDSKyXdKmle0qGpbe8KaRm+L+lFifXPknRz2PZRSVnhhANH0nsl/VTSjWE6JrEt89qqTq90GsOCpLvCPXRjHM4qaYOkyyT9IHzuPuh+dkLSmZLul3RLYl3H/g/r/bZSceFSDLcALwe+nlwp6WCi6I6nAkcBnwzpHCBK0XASC2kbigxfLJp/NLNDwnQx9Ly2ypJIp3E0cDBwfLiWYeXI8LvELzXvBK4ws4OAK8JyVTmb9vs+s//Der+tZFy4FICZ3W5m38/YdCxwvpk9ZmY/IormOEzS3sB6M/uf4Iw7F3hpH7tcBJnXNuA+5aGZTsPMZoA4ncaocCxwTpg/hwrfV2b2deDB1OpO/R/W+23F4sKlXDqlZdgnzKfXV5U3h2yrZybMFJ2ureoMa7+zMOBrkq6XdFJYt1c8xiF8Pm5gvVsanfo/Sr/bisBDkXMi6XLg8Rmb/srMvtTpsIx11mX9QOh2bUTmu/cT9e/9wAeBN1Kxa1gEw9rvLJ5rZvdIehxwmaTvDbpDJTJKv9uKwIVLTszshUs4rFNahi1hPr1+IOS9NkmfAb4SFiuZciIHw9rvNszsnvB5v6QLicxE98UZcYP5dQBpgZdFp/6PzO+2UnCzWLlcBBwnaZWkA4gc998O6v4jkg4PUWInAJ20n4ES/uAxLyMKXoAO19bv/i2BPOk0Ko+ktZJ2ieeB3yb6bS4CTgy7nUhF76sudOr/sN5vKxbXXApA0suAjwF7Al+VdKOZvSikabiAqPbCLHCymcXFtd9EFC2zBvjPMFWR0yUdQmSCuAv4I4Ae11ZZOqXTGHC3lsJewIUhgn0C+KyZXSLpOuACSb8P/Bh45QD72BVJnwOOADZK2gKcCvwdGf0f1vttJeMj9B3HcZzCcbOY4ziOUzguXBzHcZzCceHiOI7jFI4LF8dxHKdwXLg4juM4hePCxXEcxykcFy5ORySdJmnRmQkk7S/pNYnlQyV9tNjeVQtJU+HzCZK+sIzznCJpcpHHXBXS0L8ksXxor+M6nOvXJd2WTIPvOEvBhcsKp1vacjN7j5ldvoTT7g80hYuZbTaztyzhPJVEUsfBx2Z2j5m9YhmnPwVYlHAJvNbMlp1pwMz+m6geu+MsCxcuI0rQHr4n6ZyQ0fgL8RtxKDL1HknXAK+UdIika8N+F8aZjyWdLekVYf5ZkmzmSEcAAAR2SURBVK4OGXgvjdPCSPpFSZdL+q6kGyQdSDTK+tdDEas/k3SEpK+E/TdI+o/Q1rWSnh7WvzdkXb5K0p2SOgojSb8ajl8d0qDcKulpXfZ/h6KiWt+V9HdhXadr7rT+Kkl/K+lq4K0hfcz/SLpO0vtT3/stYf71kr4o6RJFxa9OT+z3KUmbQ9/fF9a9BXgCcKWkK8O63w7t3CDp85LWLeIeGAu//9+E5SlJfx9+w8slHZb4vl+S97yOkwsz82kEJyLtwYgy5wKcCbw9zN8FvCOx703Ab4b504APh/mzgVcANeCbwJ5h/auJ0qYAfAt4WZhfTfTWfQTwlcT5m8tEaXJODfPPB24M8+8NbawCNgIPALUu1/c3wD8QFf56V5f9jg7nnQzLG3pcc6f1VwGfTJz3IuCEMH8yMJX43m8J868H7gR2Dd/N3cB+qX6Mh3M/PfHbbAzzG4kK0K0Ny/8f8J6Ma7wKODS1fDjwOaKs3fF6A44O8xcCXwu/7TPi3yF9DT75tNTJNZfR5idm9o0w/6/A8xLb/g1A0q7AbmZ2dVh/DvAbqfM8BXgaUVr3G4F3A/sqSpy4j5ldCGBmj5rZdI8+PQ/4l7D/fwF7hD4AfNWiYlBbibLh7tXlPKcBvwUcCpzeZb8XAmfF/TKzBztdc47v4t8S888lengTX08HrjCzbWb2KFFerCeF9a+SdAPwHaLqilnVMA8P678RvvcTE8f34tNEAuIDiXUzwCVh/mbgajNrhPn9c57XcXLhiStHm3TiuOTyjkWcR8CtZvaclpXS+iX0qVtdjscS6+bofn9uANYRvXmvpvP1iOLqfqTbyHPetmtSlNX37cCvmtlDks4muoY0Ai4zs+OX0NdvAkdK+mAQbAANM4v7PB/3zczmu/mRHGcpuOYy2jxRUiwQjgeuSe9gZtuAhyT9elj1OuDq1G7fB/aMzyWpJumpZrYd2CLppWH9quDXeQTYpUOfvg68Nux/BLA1nGexbAL+GjgP+Psu+30NeGPC37Sh0zXn/C5ivkGUrp/4ehbBeiJBtU3SXkSmu5jkd3ct8FxJvxj6PinpyTnb+GfgYuDzLjicQeA33WhzO3CipE8DPyCqKpnFicAZ4QF8J/CGxDYzs5ng2P9oMB1NAB8GbiV6AH9a0mlAgyhF+k3ArKTvEvltvpM433uBsyTdBEyzULsjN5JOAGbN7LOKot2+Ken5wczWgkVp6A8BNkuaIXrg/mWXa+72XSR5K/BZSW8F/n0x/Tez70r6DtH3dyeRoIrZBPynpHvN7EhJrwc+J2lV2P5u4H9ztvOh8Hv9i6TFCkDHWRaecn9EkbQ/kRO9YxRVjnN8GfiQmV1ZVL+c4pF0FVGwxuaCzrc/y7x3HMfNYk4mks4kivxqM6U5leNB4OwiwomDSfDLwNZl98pZ0bjm4lQWSXsAV2RseoGZPZDa95dpj9p6zMyeXVb/HMfpjAsXx3Ecp3DcLOY4juMUjgsXx3Ecp3BcuDiO4ziF48LFcRzHKRwXLo7jOE7h/D8CMGOBCDcVwgAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x288 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"bom.precipitation.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Check Projection"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<xarray.DataArray 'proj' ()>\n",
"array(-127, dtype=int8)\n",
"Attributes:\n",
" grid_mapping_name: albers_conical_equal_area\n",
" semi_major_axis: 6378137.0\n",
" inverse_flattening: 298.257222101\n",
" standard_parallel: [-18. -36.]\n",
" longitude_of_central_meridian: 144.75200000000004\n",
" latitude_of_projection_origin: -37.852000000000004\n",
" false_easting: 0.0\n",
" false_northing: 0.0"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bom.proj"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### feed proj to cartopy\n",
"\n",
"Since the data is supplied in some azimuthal equidistant projection (or the like), we can use cartopy to do the transformation for us on the fly."
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"proj_attrs = bom.proj.attrs.copy()\n",
"proj_attrs['central_longitude'] = proj_attrs.pop('longitude_of_central_meridian')\n",
"proj_attrs['central_latitude'] = proj_attrs.pop('latitude_of_projection_origin')\n",
"proj_attrs['standard_parallels'] = proj_attrs.pop('standard_parallel')\n",
"proj_attrs.pop('grid_mapping_name')\n",
"globe_kwargs = {'semimajor_axis': proj_attrs.pop('semi_major_axis'),\n",
" 'inverse_flattening': proj_attrs.pop('inverse_flattening') }\n",
"proj_attrs['globe'] = ccrs.Globe(**globe_kwargs)\n",
"\n",
"aeqd = {}\n",
"aeqd['central_longitude'] = bom.proj.attrs.get('longitude_of_central_meridian')\n",
"aeqd['central_latitude'] = bom.proj.attrs.get('latitude_of_projection_origin')\n",
"\n",
"albers = ccrs.AlbersEqualArea(**proj_attrs)\n",
"aeqd = ccrs.AzimuthalEquidistant(**aeqd)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### get pysteps colormap, norm and clevels"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"cmap, norm, clevs, clevsStr = pysteps.visualization.precipfields.get_colormap('depth', 'mm', 'pysteps')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### xarray/cartopy powered plotting"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<cartopy.mpl.gridliner.Gridliner at 0x7f71d3e5dbe0>"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAnkAAAIDCAYAAAB4sAgGAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzs3Xu8HWV97/HvLzt3IAkhBmMAuUViNBK2XGLFCnLoCYSK9ULBY0UrpahYbWtrtG1sYz0n3trSCtIYUWytNFqq0eSF9UWNSmtiwm7KFggSkMsmgRBoAiEx19/5Y63ZmT2ZWdeZWWvNfN77tV97z6yZeZ41a+1Z3/08M8+YuwsAAADFMqrTFQAAAED6CHkAAAAFRMgDAAAoIEIeAABAARHyAAAACmh0rQdnzJjhTz75ZF51AQAAnfGou5/c6UogXVZrCBUz84MHD2rUKBr88jY4OKgZM2Zo2rRpna5KKa1Zs0YXXHBBp6tRSuz7zmHfd866des0d+5cTZw4sSPlm5nc3TpSODJDegMAACggQh4AAEABEfIAAAAKiJAHAABQQIQ8AACAAiLkAQAAFBAhDwAAoIAIeQAAAAVEyAMAACggQh4AAEABEfIAAAAKaHS9Bfbu3cu9azvg4MGD2rdvn/bu3dvpqpTSoUOH2Pcdwr7vHPZ95wTH/L6+vk5XBQVi7p78oJn/x3/8h8y4Z3Henn/+eY0bN05jx47tdFVK6dlnn9XUqVM7XY1SYt93Dvu+c3bs2KFjjjmmYyHvV37lV+TufNgXTN2Qd/DgQVryOmBwcFAzZszQtGnTOl2VUlqzZo0uuOCCTlejlNj3ncO+75x169Zp7ty5mjhxYkfKNzNCXgGR3gAAQFcys3PN7IJO16NX1T0nDwAAIG9mZscff/w/SZpgZie7+/5O16nX0JIHAAC6zvjx49+0YMGC466++upJkyZNuqbT9elFhDwAANBVzMwmT578mU984hNTPvKRjxw9YcKEPzWzMZ2uV68h5AEAgK4StOKdeOKJmjp1qmjNaw0hDwAAdI1wK14wj9a81hDyAABA1wi34gVozWtNx8bJ+9LHb0x9m0Vy1MzJ2rtzjw7s2tfpqpTS5DOma+cD2zpdjVJi33fOlNnTtWMT+74Tjjn1OL3w+A4d2n8wcZn3/MX7Myu/W8bJq15R++D69etPC4c8qTJY95w5c7Y89dRTXGnbIFryAABAV4hrxQvQmtc8Qh4AAOi4uHPxojg3rzmEPAAA0HG1WvECtOY1h5AHAAA6qpFWvACteY0j5AEAgI5qpBUvQGte4wh5AACgY5ppxQvQmtcYQh4AAOiYZlrxArTmNYaQBwAAOqKVVrwArXn1EfIAAEBHtNKKF6A1rz5CHgAAyF07rXgBWvNqI+QBAIDctdOKF6A1rzZCHgAAyFUarXgBWvOSEfIAAECu0mjFCzTTmmdmt5jZNjP7WWT+B8zsATO718w+HZr/UTPbXH3sfydsc6qZfd/MHqz+PLbtJ5USQh4AAMhNmq14gSZa874iaUGkPhdKulzSq9z9FZI+W50/R9KVkl5RXecmM+uL2eYiSXe6+yxJd1anuwIhDwAA5OnX02rFC0ydOlXvete7Jkt6Z63l3P1Hkp6NzH6vpKXuvre6zLbq/Msl3ebue939F5I2Szo3ZrOXS7q1+vutkt7U2rNI3+hOVwAAAJTKeReuvHXKYy+5tf6STTh/n476lHS1mb03NHuZuy+rs+rLJL3OzD4p6ZeSPuzu6yXNlLQ2tNxQdV7U8e6+VZLcfauZTW/9WaSLkAcAAHJlWWyzstFN7n5tk6uOlnSspPmSzpG0wsxOVXw1vZ065o3uWgAAUGZDkm73ip9KOiRpWnV+uE/5BElbYtZ/ysxmSFL157aYZTqCkAcAAMrsW5LeIElm9jJJYyVtl7RS0pVmNs7MTpE0S9JPY9ZfKenq6u9XS/p25jVuECEPAACUgpl9XdJPJJ1hZkNm9h5Jt0g6tTqsym2Srq626t0raYWk+yTdIen97n6wup3lZnZ2dbNLJV1sZg9Kurg63RU4Jw8AAJSCu1+V8NA7Epb/pKRPxsy/JvT7M5IuSqWCKaMlDwAAoIAIeQAAAAVEyAMAACggQh4AAEABEfIAAAAKiJAHAABQQIQ8AACAAiLkAQAAFBAhDwAAoIAIeQAAAAVEyAMAACggQh4AAEABEfIAAAAKiJAHAABQQIQ8AACAAiLkAQAAFBAhDwAAoIAIeQAAAAVEyAMAACggQh4AACgFM7vFzLaZ2c9C86aa2ffN7MHqz2NDj33UzDab2QNm9r8Ttpm4fqcR8gAAQFl8RdKCyLxFku5091mS7qxOy8zmSLpS0iuq69xkZn0x24xdvxsQ8gAAQCm4+48kPRuZfbmkW6u/3yrpTaH5t7n7Xnf/haTNks6N2WzS+h1HyAMAAEUx28w2hL6vbWCd4919qyRVf06vzp8p6fHQckPVeY2u33GjO10BAACAlGxy90aCXSMsZp6ntO1c0JIHAADK7CkzmyFJ1Z/bqvOHJJ0YWu4ESVuaWL/jCHkAAKDMVkq6uvr71ZK+HZp/pZmNM7NTJM2S9NMm1u84Qh4AACgFM/u6pJ9IOsPMhszsPZKWSrrYzB6UdHF1Wu5+r6QVku6TdIek97v7wep2lpvZ2dXNxq7fDTgnDwAAlIK7X5Xw0EUJy39S0idj5l8T+v2ZpPU7rW7I+/GPf5xJwZPP6JqLT7pS39g+jZk0Xn6op87xLIy+8aN5j3YI+75z+sax7zulb9xoTTrtOHmNQ/4Pf/jD/CqEQqgb8l772tdq1Kj0e3W//BdfSH2bRXLUzMnau3OPDuza1+mqlNLkM6Zr5wNdc+5sqbDvO2fKbPZ9pxxz6nF64fEdOrT/YOIyb77yLTnWCEVQN+SNGjUqk5AHAAAax2cxmsU7BgAAoIAIeQAAAAVEyAMAACggQh4AAEABEfIAAAAKiJAHAABQQIQ8AACAAiLkAQAAFBAhDwAAoIAIeQAAAAVEyAMAACggQh4AAEABEfIAAAAKiJAHAABQQIQ8AACAAiLkAQAAFBAhDwAAlIaZTTGzb5rZJjO738xeY2afMLN7zGyjmf2bmb0kYd0FZvaAmW02s0V5171ZhDwAAFAmN0i6w91nSzpT0v2SPuPur3L3eZK+K2lxdCUz65N0o6RLJM2RdJWZzcmv2s0j5AEAgFIws0mSflXSlyTJ3fe5+w53fy602FGSPGb1cyVtdveH3X2fpNskXZ51ndsxutMVAAAA5TJ9gXTS6elu8+mtkpZrtpltCM1e5u7LQtOnSnpa0pfN7ExJd0v6oLu/YGaflPROSTslXRhTxExJj4emhySdl+ZzSBsteQAAoCg2ufvZoe9lkcdHS+qX9AV3P0vSC5IWSZK7/4m7nyjpa5Kuj9m2xcyLa/HrGoQ8AABQFkOShtx9XXX6m6qEvrB/kvSWhHVPDE2fIGlL6jVMESEPAACUgrs/KelxMzujOusiSfeZ2azQYm+UtClm9fWSZpnZKWY2VtKVklZmWuE2cU4eAAAokw9I+lo1qD0s6d2SlleD3yFJj0q6TpKqQ6ksd/dL3f2AmV0v6XuS+iTd4u73duQZNIiQBwAASsPdN0o6OzI7rntW7r5F0qWh6dWSVmdXu3TRXQsAAFBAhDwAAIACIuQBAAAUECEPAACggAh5AAAABUTIAwAAKCBCHgAAQAER8gAAAAqIkAcAAFBAhDwAAIACIuQBAAAUECEPAACggAh5AAAABUTIAwAAKCBCHgAAQAER8gAAAAqIkAcAAFBAhDwAAIACIuQBAAAUECEPAACUipn1mdl/mdl3q9NTzez7ZvZg9eexCestMLMHzGyzmS3Kt9bNI+QBAICy+aCk+0PTiyTd6e6zJN1ZnR7BzPok3SjpEklzJF1lZnNyqGvLCHkAAKA0zOwESQslLQ/NvlzSrdXfb5X0pphVz5W02d0fdvd9km6rrte1Rne6AgAAoGRmSDot5W32SZJmm9mG0Nxl7r4ssuTfSPpjSceE5h3v7lslyd23mtn0mBJmSno8ND0k6bx2q50lQh4AACiKTe5+bdKDZnaZpG3ufreZXdDkti1mnje5jVwR8gAAQFm8VtIbzexSSeMlTTKzf5T0lJnNqLbizZC0LWbdIUknhqZPkLQl8xq3gXPyAABAKbj7R939BHc/WdKVkv7d3d8haaWkq6uLXS3p2zGrr5c0y8xOMbOx1fVX5lDtlhHyAABA2S2VdLGZPSjp4uq0zOwlZrZaktz9gKTrJX1PlStzV7j7vR2qb0PorgUAAKXj7mskran+/oyki2KW2SLp0tD0akmr86lh+2jJAwAAKCBCHgAAQAER8gAAAAqo7jl5P//5z2UWNzRMeya8+Jj6C5XY6KPGSn2mMUeP63RVSmnUmFG8RzuEfd85Npp93yl940Zr/PSj5QcPJS7zwAMP5FgjFEHdkHfMMcdo1Kj0G/wO7tmf+jaL5NDEsTq094AO/vJAp6tSSn7IeY92CPu+g5x93yk++ZAO/nK//EByyJs0aVKONUIR1A15M2bMyCTk7dv5y9S3WSRjjh6n/S/s04Fd+zpdlVKa8OJJvEc7hH3fORNnsO87ZdxxR2n/c3t1aP/BxGVmzJiRY41QBJyTBwAAUECEPAAAgAIi5AEAABQQIQ8AAKCA6oY8d8+jHgAAAEgRLXkAAAAFRMgDAAAoIEIeAABAARHyAAAACoiQBwAAUECEPAAAgAIi5AEAABQQIQ8AAJSCmZ1oZj8ws/vN7F4z+2B1/p+b2RNmtrH6fWnC+gvM7AEz22xmi/KtffNGd7oCAAAAOTkg6Q/dfcDMjpF0t5l9v/rYX7v7Z5NWNLM+STdKuljSkKT1ZrbS3e/LvNYtoiUPAACUgrtvdfeB6u/PS7pf0swGVz9X0mZ3f9jd90m6TdLl2dQ0HYQ8AABQOmZ2sqSzJK2rzrrezO4xs1vM7NiYVWZKejw0PaTGA2JH0F0LAADydaKk2ZlsebaZbQhNL3P3ZdGFzOxoSf8i6UPu/pyZfUHSJyR59efnJP12dLWY8jydameDkAcAAIpik7tfW2sBMxujSsD7mrvfLknu/lTo8S9K+m7MqkOqxNPACZK2tF3jDNFdCwAASsHMTNKXJN3v7n8Vmj8jtNhvSPpZzOrrJc0ys1PMbKykKyWtzLK+7aIlDwAAlMVrJf2WpEEz21id9zFJV5nZPFW6Xx+R9LuSZGYvkbTc3S919wNmdr2k70nqk3SLu9+b9xNoBiEPAACUgrvfpfhz61YnLL9F0qWh6dVJy3YjumsBAAAKiJAHAABQQIQ8AACAAiLkAQAAFBAhDwAAoIAIeQAAAAVEyAMAACggQh4AAEABEfIAAAAKiJAHAABQQIQ8AACAAiLkAQAAFBAhDwAAoIAIeQAAAAVEyAMAACggQh4AAEABEfIAAAAKiJAHAABQQIQ8AACAAiLkAQCAUjCzW8xsm5n9LDL/A2b2gJnda2afTlh3QXWZzWa2KJ8at4eQBwAAyuIrkhaEZ5jZhZIul/Qqd3+FpM9GVzKzPkk3SrpE0hxJV5nZnMxr2yZCHgAAKAV3/5GkZyOz3ytpqbvvrS6zLWbVcyVtdveH3X2fpNtUCYZdbXSnKwAAANBBL5P0OjP7pKRfSvqwu6+PLDNT0uOh6SFJ57VSmJn1x8zeKelRdz/QyjaTEPIAAEC+ZqvFiFRDJdHMNrMNobnL3H1ZA2seK2m+pHMkrTCzU93dQ8tYzHoeM68RN0nql3RPdbuvrP5+nJld5+7/1uJ2j0B3LQAAKIpN7n526LtewJMqrXK3e8VPJR2SNC1mmRND0ydI2tJiHR+RdFa1fq+WdJakn0n6X5JiL/poFSEPAEpu7ZKNWrtko3Zt3a21SzZ2ujpA3r4l6Q2SZGYvkzRW0vbIMuslzTKzU8xsrKQrJa1ssbzZ7n5vMOHu96kS+h5ucXuJCHkAUGKEOpSJmX1d0k8knWFmQ2b2Hkm3SDq1OqzKbZKudnc3s5eY2WpJqp4rd72k70m6X9KKcFBr0gNm9gUze331+yZJPzezcZL2t/kUR+CcPAAoKQIeysbdr0p46B0xy26RdGloerWk1SlU412S3ifpQ6qck3eXpA+rEvAuTGH7wwh5AFBCa5ds1CIdHs91qZZ2sDZAqcxx989J+lwww8x+3d2/I2lXmgXRXQsAJTTQf42u6N+u0685TVf0Hz79aP7ieR2sFVAKXzSzucGEmV0p6U+zKIiQB3S5C5Ys0E3fPVvHLZnZ6aqgQN532QatGJim/oF1WjFQuZDw6BkTO1wroBTeKulWM3u5mf2OpPdL+rUsCiLkAV3sgiULhltZVmlVh2uDolmz+A71DyzXUi3t6hY8zh1EkVSvor1S0r+oEvh+zd13ZlEW5+QBXWyplmrFQOW8qTVd/CGM3tXN4U46HPDWLtnY9XUFajGzQY0cQHmqpD5J68xM7v6qtMsk5AFdbP7ieVqjOzpdDSAzQYiLC3C04KFgLsu7QLprAQAdE4Q7Ah2Kzt0frfWdRZl1W/J27typUaPSz4J9E8akvs0isdGj1DdutPxgq7fGQ1tGGe/RTmHfd47lv+8Hlz2gSadMknTk50IwP1Dk94X1mfrGj5aNTv683bFjR441QtrMbMDd+9tdphl1Q94jjzwis7j78rZnwouOTn2bRTJ64hjZ6FEac9S4TlellPpGj+I92iHs+84ZNSb/fX/ypScP//780F5NP+s4SdK2/3pmxGPRx4umb+xojTvuKOlQ8j/2v/jFL3KsETLwcjO7p8bjJmlymgXWDXlnnnmm+vr60ixTkjTwrf9MfZtFctTMydq7c48O7NrX6aqU0uQzpmvXY//T6WqUEvu+c6ZMzH/f33PjPVqohZqruRrUoFZpleYvnqeJxx1u0Qp35Rb14otjTj1Ou5/YqUP7DyYuc9Z7zsqxRsjA7AaWSX4DtIALL4ACuObFo6T+5Vq++rc7XRWgYeHwFgS86Pz5i+dp/uJ5nLOHnpfVeXe1cOEF0OOuufQWXfPe90kD13S6KkBTgla5VdWvWq10RW3BA7JEyAN62NolGw+Hu4H3drYyQAuClrogxNULc7ToAY0j5AG9rn+5JOlfl3+rwxUBshUeboWwB9RHyAN62CJV7oax/As36ZnFT3S4NkD6omEu2vIH9Boze7OZPWhmO83sOTN73syey6IsQh7QwwY1qP733aTlTx7qdFUAAI35tKQ3uvtkd5/k7se4+6S6a7WAq2uBHvbM4if0PtGCBwA95Cl3vz+PgmjJAwDkqpnz6eiWRQFtMLN/NrOrql23bzazN2dREC15AHJzwZIFOv2a04an+wfWacXANK1ZfEcHa4WsrV2ycURYSwpuXEyBkpgkabekXwvNc0m3p10QIQ9AbtYsvkPXLNmoRVqkK/q3a6D/PJ3eL10TCQEollZeW94PyIKZ3SLpMknb3P2V1XlTJf2zpJMlPSLpCnc/4tYvZrZA0g2S+iQtd/elrdTB3d/dUuVbQMgDkKv5i+dpje7Q+yRJGyoz+UCHCHbIxVckfV7SV0PzFkm6092Xmtmi6vRHwiuZWZ+kGyVdLGlI0nozW+nu9zVbATM7QdLfSXqtKi14d0n6oLsPNf90auOcPAAAUAru/iNJz0ZmXy7p1urvt0p6U8yq50ra7O4Pu/s+SbdV12vFlyWtlPQSSTMlfac6L3WEPAAAUBSzzWxD6PvaBtY53t23SlL15/SYZWZKejw0PVSd14oXufuX3f1A9fsrkl7U4rZqorsWAADkapOkYzPYpqRN7t5IsGuWxczzFre13czeIenr1emrJD3T4rZqoiUPAACU2VNmNkOSqj+3xSwzJOnE0PQJkra0WN5vS7pC0pOStkp6a3Ve6mjJAwAAZbZS0tWSllZ/fjtmmfWSZpnZKZKekHSlpLe3Upi7Pybpja1VtTmEPAAAUApm9nVJF0iaZmZDkj6uSrhbYWbvkfSYpLdVl32JKkOlXOruB8zseknfU2UIlVvc/d4my/5jd/+0mf2dYrp63f332nhqsQh5AACgFNz9qoSHLopZdoukS0PTqyWtbqP44FZmG9rYRlMIeQAAABlz9+9Uf93t7t8IP2Zmb8uiTC68AAAAyM9HG5zXNlryAAAAMmZml6jS/TvTzP429NAkSQeyKJOQBwAAkL0tqpyP90ZJd4fmPy/p97MokJAHAACQMXf/b0n/bWb/5O778yiTkAcAAJCfk83s/0maI2l8MNPdT027IC68AAAAyM+XJX1BlfPwLpT0VUn/kEVBhDwAAID8THD3OyWZuz/q7n8u6Q1ZFER3LQAAQH5+aWajJD1YvYvGE5KmZ1EQLXkAAAD5+ZCkiZJ+T9KrJb1D0juzKIiQBwAAkJ+T3X2Xuw+5+7vd/S2STsqiILprAaCHrF2ycfj3+YvndbAmAFr0UUnfaGBe2wh5ANDlwsEuOp+gB/QG7ngBAABQTNzxAgCQbKEWapVWDU8ft2SmVmkVLXpAlwvd8eJr7p5Jy10UIQ8Autz8xfNGdNku1EJJ0qrql8S5ekC3M7MV7n6FpP8yM48+7u6vSrtMQh4A9IhoK154vqQRgY+gB3SdD1Z/XpZXgQyhAgA9IingBfMXVr+kwy17a5dsTLxwA0B+3H1r9eejkvZKOlPSqyTtrc5LHSEPAHrA/MXzYlvnVmnVcLALRIMegO5hZtdI+qmkN0t6q6S1ZvbbWZRFyAOAHhINe3Hn58Uh8AEVZjbFzL5pZpvM7H4ze03k8cvN7B4z22hmG8zs/JSr8EeSznL3d7n71arc9eIjKZchiXPyAKAnBRdjRINdUtDjHD1g2A2S7nD3t5rZWFVuMRZ2p6SV7u5m9ipJKyTNTrH8IVWGTQk8L+nxFLc/jJY8AOhRBDegOWY2SdKvSvqSJLn7PnffEV6mesux4OrXoyQdcSVsm56QtM7M/tzMPi5praTNZvYHZvYHaRZEyAOAHkbQA0aYXe1iDb6vjTx+qqSnJX3ZzP7LzJab2VHRjZjZb5jZJkmrJKV9vtxDkr6lw+Hx25K2Sjqm+p0aumsBoMeFg1703DtCILrRZknjUt7mY5Ufm9w9GuzCRkvql/QBd19nZjdIWiTpz8ILufu/SvpXM/tVSZ+Q9L/Sqqe7/0Va26qHkAcAAMpiSNKQu6+rTn9TlZAXy91/ZGanmdk0d9/eTsFm9jfu/iEz+45iuoDd/Y3tbD8OIQ8ACoJWPKA2d3/SzB43szPc/QFJF0m6L7yMmZ0u6aHqhRf9ksZKeiaF4v+h+vOzKWyrIYQ8ACgQgh1Q1wckfa16Ze3Dkt5tZtdJkrvfLOktkt5pZvsl7ZH0m6ELMVrm7ndXf90gaY+7H5IkM+tT+r3Xkgh5ANCzoverJeAB9bn7RklnR2bfHHr8U5I+lWEV7lTlHL9d1ekJkv5N0q+kXRAhDwB6CIMaAz1vvLsHAU/uvsvMomP1pYKQBwA9gnPugEJ4wcz63X1Akszs1ap0C6eOkAcAPYCABxTGhyR9w8y2VKdnSPrNLAoi5AFADyDUAcXg7uvNbLakMySZKmP77Y9b1sxWNrDJZ939XXEPEPIAAAByUj3/7g8kvdTdf8fMZlWHdPluzOIvl3RNrc1JujHpQUIeAABAfr4s6W5Jr6lOD0n6hqS4kPcn7v7DWhszs8Q7aHDvWgAAgPyc5u6flrRfktx9jyotckdw9xX1NlZrGUIeAABAfvaZ2QRVb21mZqdJ2hu3oJn1mdnvmtknzOy1kcf+tF5BhDwAAID8fFzSHZJONLOvqTI48h8nLPv3kl6vym3V/tbM/ir02JvrFcQ5eQAAADkwM5O0SZWANl+VbtoPuvv2hFXOdfdXVdf9vKSbzOx2SVcpoYs3jJY8AACAHFTvgfstd3/G3Ve5+3drBDxJGhta94C7Xytpo6R/l3R0vfLqtuT9+Mc/bqDazZt8xvRMtpsFq5uV0zdqTJ/GTB4vHWr7nshoQd+40Zoyu3feo0XCvu+cvvHs+2a1f9v6ir5xozXptONqbm/NmjXpFIZOW2tm57j7+gaW3WBmC9z9jmCGuy+pDqT8hXor1w15r3vd69TX19dAPZrzpY8nDusCSUfNnKy9O/fowK59na5KKU0+Y7p2PrCt09UopbLv+7VLNnZs4OMps6drx6by7vtOOubU4/TC4zt0aP/BxGXe+va35VgjZOhCSdeZ2SOSXlCl29WDbtkwd39H3Abcfbmk5fUK4pw8AOgSnQx4AHJzSTsrm9myardtXYQ8AOgC0XvTAigmd3/UzPolna/KMCr/4e4DTWzinEYXJOQBQAeFwx2teEDxmdliSW+TdHt11pfN7Bvu/pcNbuLJRssi5AFAFyDgAaVxlaSz3P2XkmRmSyUNSEoMeWZ2tqQ/kfRSSWPMbFAJ5/GFEfJQaBcsWaA1i++ovyDQIYQ7oHQekTRe0i+r0+MkPVRnna9J+iNJg5IONVoQIQ+FtXbJRt3Uf7be1+mKAABw2F5J95rZ91U5J+9iSXeZ2d9Kkrv/Xsw6T7v7ymYLIuShsBZqoZb3n6e1S5bTWgIAGGZmfZI2SHrC3S+LPGaSbpB0qaTdkt7V5IUR9fxr9TuwpoF1Pm5my1W5BdrwfW7d/fbkVQh5KIGFWqhn9ESnqwEA6B4flHS/pEkxj10iaVb1+zxVBh0+L62C3f3WFlZ7t6TZksbocHet6/DFG7EIeSisVVql/oEZna4GAKCLmNkJkhZK+qSkP4hZ5HJJX63egmytmU0xsxnuvjXPekac6e5zm12Je9ei8Oaq6b8LAEBvmm1mG0LfcYMG/42kP1byBQwzJT0emh6qzuuktWY2p9mVaMkDAAC52qzQiWUp2V75sanW3SDM7DJJ29z9bjO7IGmxmHmdvpH8+ZKuNrNfqLLrEm+FFkbIQ+Fd0b+dK2wBAJL0WklvNLNLVRnGZJKZ/WPkHrFDkk4MTZ8gaUu7BZvZd1QjLLr7G2usvqCVMgl5AACgFNz9o5I+KknVlrwPRwKeJK2UdL2Z3abKBRc7Uzof77PVn2+W9GJJ/1idvkqVsfNq1fvRVgok5KGw5i+ep4HvdroWAIBuZ2bXSZK73yxptSrDp2y6tB3iAAAgAElEQVRWZQiVd6dRhrv/sFrWJ9z9V0MPfcfMfpRGGVGEPBTaioFpuqJ/e6erAQDoMu6+RtUx6qrhLpjvkt6fYdEvMrNT3f1hSTKzUyS9KIuCuLoWhbZUS7ViYFqnqwEAQOD3Ja0xszVmtkbSDyR9KIuCaMlDoc1fPE9rxL1rAQDdwd3vMLNZqgxuLFWuCE77YmNJtOQBAADkxswmSvojSde7+39LOqk6tEvqCHkAAAD5+bKkfZJeU50ekvSXWRREyAMAAMjPae7+aUn7Jcnd9yh+AOa2EfIAAADys8/MJqg6MLKZnab0bwAiiQsvAAAA8vTnku6QdKKZfU2Vu3CkMhZfFCEPAAAgJ+7+b2Z2t6T5qnTTftDdMxnQle5aAACAnJjZne7+jLuvcvfvuvt2M7szi7JoyQMAAMiYmY2XNFHSNDM7Vocvtpgk6SVZlEnIAwAAyN7vqnJni5dIuluHQ95zkm7MokBCHgCgdNYu2Tj8+/zF8zpYE5SFu98g6QYz+4C7/10eZRLyAAClFg58EqEP2XL3vzOzV0qaI2l8aP5X0y6LkAcAAJATM/u4pAtUCXmrJV0i6S5JqYc8rq4FAJRKtOUOyNlbJV0k6Ul3f7ekMyWNy6IgQh4AACGEQGRsj7sfknTAzCZJ2ibp1CwKorsWAIAQzslDxjaY2RRJX1TlKttdkn6aRUGEPABAadBKh05z9/dVf73ZzO6QNMnd78miLEIeAKA0glY6wh7yZmb9tR5z94G0yyTkAQBKJ9olS+grBzM7UZWrWF8s6ZCkZdXx68LLXCDp25J+UZ11u7svSaH4z9V4zCW9IYUyRiDkAQBKj/PwSuOApD909wEzO0bS3Wb2fXe/L7Lcj939sjQLdvcL09xeIwh5AACgFNx9q6St1d+fN7P7Jc2UFA15mTGzdybUjcGQAQAAEsw2sw2h6WXuvixuQTM7WdJZktbFPPwaM/tvSVskfdjd702xjueEfh+vyph5A8pgMGRCHgAAyNX3Np0mHTs53Y1u2i1p0yZ3v7beomZ2tKR/kfQhd38u8vCApJe6+y4zu1TStyTNSqua7v6BSF0mS/qHtLYfxmDIAACgNMxsjCoB72vufnv0cXd/zt13VX9fLWmMmU3LsEq7lWKIDCPkldjaJRuHryjjyjIAQNGZmUn6kqT73f2vEpZ5cXU5mdm5qmSlZ1Ksw3fMbGX1+7uSHlDlat7U0V1bUuFQR8ADAJTEayX9lqRBMws+/D4m6SRJcvebVbm37HvN7ICkPZKudHdPsQ6fDf1+QNKj7j6U4vaHEfJKKCnUrV2ykWEEAACF5e53SbI6y3xe0uczrMMPJal639rR1d+nuvuzaZdFdy1GoFUPAIDsmNm1ZvaUpHskbVDl/rUbaq/VGlrySiYIcfMXzxtuuSPYAUA2osdXeksg6Y8kvcLdt2ddECGv4BrpmuWgAwDpCx9/Oc4i5CFVrqjNHCGvgJKC3UIt1CqtOmI5Dj4AkD6OrUjwUUn/aWbrJO0NZrr776VdECGvQILQtlALJUlzNfeIZcLzBjWoVVpF2AMAID9/L+nfJQ1KOpRlQYS8gggHvGi4u6L/yG7/FQPTNLf6tVRLM6mLRHAEACDigLv/QR4FEfJ6XDhQhQNeXLALu6J/u1YMTBteL9yNm0ZdAABArB+Y2bWSvqOR3bWpD6FCyOtRcYEqLuAN9J+n/oF1Gug/b8Sy/QPrRgS9LOoDAACO8Pbqz4+G5rmkU9MuiJDXg6KBKjgHTzoy4IV/NrLdZrtXoy2JgTRaBgEAKBp3PyWvsgh5PSYu4CW14OVZl3DAiy7DeXkAAFSY2Tvj5rv7V9Mui5CXsTSvXK0V8MKOCHj9Xwg9+N626xGtS1LAAwAARzgn9Pt4SRdJGpCUesjjtmYZCgehtUs2tnXeWq0uWulwK17NgBeSZUtftKuW8/UAAKhw9w+Evn9H0lmSxmZRFi15KWk0yDTTsnffrZv1sreclBjwot20scFt4L2xQa9/YN3w78F2gjHzmm11rHf+Hd21AAAk2i1pVhYbJuS1qN3WqVr3M1y7ZKPmvHtO7HJJAa8hMV21wRW2czW3qYslGr3nLQEPAIDDzOw7qlxNK1V6VOdIWpFJWe6e/KCZDwwMaNSo9Ht171//s9S3mbVtdz+TW1nHnnGsdj+9W3ufHR5CRyfpJEnSVE0dnvf9ab8c/n37i45vePvTnn5qxPTF28dLkp7Vs3pMj2n6q49rqr5J+6bZ7XSLsVMmaN+OPZ2uRimx7ztn3JQJ2su+74ixk8Zr/6698kPJn8kvP+eVmZU/b948ubtlVkCImX1Sf33ax/SayelueNNu6V2bvuju16a74XSZ2etDkwckPeruQ1mUVbcl76STTsok5G28Y139hbrE4PKf13z8fJ2f+Nhduqt+ATF/030T+vTMvc/ouYeeG54XDXmLZu+QHp84/Pim2SfVLyu0NUmaveleSdJ/SVq6aYqmaqru0l2afPL4Jral4eUHv3h4X839nZdp7zMvNLWdbjHm6HE9W/dex77vnDHHsO87pW/8aO39nz3yAwcTlznppGaO8ehij0na6u6/lCQzm2BmJ7v7I2kXVDfkTZkyRX19fWmXqwO796e+zSzU6pKsdY/YQQ1KOjIANtolum/HPr3wxAva+dDOEWVJ0gRNkCRtnTzyYDx568gw2sjFFZO3jhn+/erJL2jFwDTt1E7dtWh9S12tL/8/h4f/6ZXXOI4f8p6ufy9j33cQ+75j/KDr4J79OrQ/OeQde+yxOdYIGfqGpF8JTR+szjsnfvHWcXVtgnpXw9YKeLXmL4z5akQQDpO2Gyd8cUWSuCDIkCgAAGRmtLvvCyaqv3N1bVaauYgiHIDqBa56jwetfUn3jn2NXnPE/EENthT0arXqBbc+kyoXYvzZQGU+AxkDAJC6p83sje6+UpLM7HJJTVxF2bjSt+RlFfAaER66JM5P9JMR06u0akS5cfedXTEwLXZ+I616cRjjDgBQFGZ2i5ltM7PEqz/N7AIz22hm95rZDzOoxnWSPmZmj5nZY5I+Iul3MyinvC15aYaXWsOYxAWusHpBL2qplmqRFtXdfnh+UL/+gXWJLXrh1rxonWjRAwAUxFckfV4Jd5cwsymSbpK0wN0fM7PpaVfA3R+SNN/MjlZllJPn0y4jUMqWvFYDXivnxV3Rv334Oy2D1a801Wvpo0UPANDr3P1Hkp6tscjbJd3u7o9Vl9+Wdh3M7P+a2RR33+Xuz5vZsWb2l2mXI5U05LWrnYCVFPZqXewwf/G8ES1p9Vr9ovULt+pFw1z/wLqGu3KDi1EIfACALjXbzDaEvpsdM+9lko41szVmdreZvTODOl7i7juCCXf/H0mXZlBOebtrOyG4u0ScuOA25+rTdWDX8AU4R9xlIu4ijEYuzAi6bePCXb3u5QBduACAlm1eII17RbrbfOwxSUs3tTkY8mhJr5Z0kaQJkn5iZmvdvfaAuc3pM7Nx7r5XqoyTJ2lcitsfVsqWvG4LJ9GAN3/xPM2Y/6KG1gm32kV/D09Hw1urF2KE0aIHACiYIUl3uPsL7r5d0o8knZlyGf8o6U4ze4+Z/bak7yvhHMF2lTLktaOZ+7vGqXduXr0AWqvbdm71K6k7uV4rXSOteM2M7QcAQI/5tqTXmdloM5so6TxJ96dZgLt/WtJfSnq5pFdI+oS7fyrNMgKlDXmdas2LBqlwIGu0TsFy4bAVF+zium0buRo3LDxo8yItGg6SAAD0GjP7uqSfSDrDzIaqrWnXmdl1kuTu90u6Q9I9kn4qabm7Jw630ip3v8PdP+zufyhpl5ndmHYZUsnPyQvCUtbdjkkBqpWAl7SdaPCKtuhFx9cLWhTj6hYepJlABwAoCne/qoFlPiPpM1nWw8zmSbpK0m9K+oWk27Mop9QhL1ArYCUFwHCwCoemsHrhrt2u31ValRjEBjUYu/1wnZNEt5cUCLvp3EYuBAGyw98X0D4ze5mkK1UJd89I+mdVxsm7MKsyS9td247onSfiJLWQBeErHMDSPHgmBchWA2VceE17jL52hId04UIQIF3B3xQBD0jFJlWu2v11dz/f3f9O0sEsC6Qlr47osCWBaBdpI12ycUGr3YNnXNisFeiC+jTSDRsOeMHzC0JqNxz0CXVAtrrh7xwokLeo0pL3AzO7Q9JtkizLAmnJa9EqrarbohVutYtruWv1ABqEm+hVrmm2sIXH9OuVgNcN9QIAII67/6u7/6ak2ZLWSPp9Sceb2RfM7NeyKJOQ14Ck8BAEvbjx6ZLOiWs34NVTqxUv/Fi9W6N1U8CLBjpa8AAAvao6Bt/X3P0ySSdI2iiFbkqfIrprGxQNN0HQCAen4PeFWhg7wHEawq149bpck7qaw+p134bDaicDHsEOAFA07v6spL+vfqeOkNeiaICKC3vBcmmpF3TavVo3qlsCHoDiCf99c6oFkA1CXhuSxtnL+oAV14oX7XptpBUvGIIlWD/YZrSruRcOwL1QRwAV/AMH5IOQl4I8Akatbtq4c+viDqLhCzWCEJfUAil1NjiF6x9+zku19IhlCXhoFuO+xdu1dXdDAazdfRf9J5TXA8gGIa+H1LpnbL2u2riAl6SbDrbRUBs937Gb6oreQKA40vA/kZ87fG/qtE//iMrrjkNAmRHyekD4Io9mh02pFQyj+OBDGfA+H2ntko3Dx4kJmqA92tPQOuxHoPsxhEoPibZq1bs9WtIVuHHBr5sO2I1eQdxNdQZ6UTjgSRoOeFm34oXNXzyPv2UgI7Tk9Zhoy10jAS8aDPM8gLejkbtyAGhNNOAFeuX4AKA+Ql4PCbprmzn/Liwu4HXbf9BBK96iyLiQ0VbLbqs30EvixvmUpPN1flPb4O8Q6G6EvBjdfMusRgNeI61g3fKcAkm3a5OyG3sQAICi4py8Bq1dsrFQV4H1UlAKd1H3Ur2BbpX2HXgAdCda8iK68aAVd0AOukri6hvXitfpu1fUE717SLd3KwO9Lu5v6ugZE0fMb+R4GCzD3yjQfWjJC+nGgJckOKCGr0wL370irBdawpLqxZV3AAC0hpa8Aglav8LntNUbZqWbEOaA7tLI7REDXIgBdB9CXlUet/LJUvRgHG696+ZuWgDdrZk7U6QR9Lr5wjf0PjO7RdJlkra5+ytjHjdJN0i6VNJuSe9y94F8a5keQp56q5u2EXGtdhwkAeSh1aBX6zjMeX9I0VckfV7SVxMev0TSrOr3eZK+UP3ZkzgnrwQ4MAJoR7P/CDezfNFGLkB3c/cfSXq2xiKXS/qqV6yVNMXMZuRTu/TRklcg0W4Vwh2ArATn/iad79tIix7hrsQ2XybtXZDuNrffLWnpbDPbEJq7zN2XNbGVmZIeD00PVedtTaGGuSt9yGv0INNLgamX6gqgu8UdI8MXdzVyF55GtgmkZJO7X9vG+hYzz9vYXkfRXQsAOEJSN2rcHWmSbqUYt369gLew+hXFP6/IyZCkE0PTJ0ja0qG6tK3ULXlFbMUD8hb+O+JvpffVOi4m3TZxUIOJLXrNtNrV6wIGcrBS0vVmdpsqF1zsdPee7KqVShzymj0xmA8vYKSkVhr+VnrXrq27Ex+rFfDSQrhD1szs65IukDTNzIYkfVzSGEly95slrVZl+JTNqgyh8u7O1DQdpQt5nAsCZIuWvd60dslGnf+582MfS+qOjRuPMwu8j5AWd7+qzuMu6f05VSdzpQt5tdBVADSGf5aKpZEu2kCaLXeNIOABrePCi6qk/1SBsgsCQPQnyiHpn95Vka8sEPCA9pSqJS/pw4mAB8RrJ+DxAd0bGnlN8+7d4L0DpKM0Ia/WgSx6AIu7SowTylE27bTY8beCWnh/APmguzYiaRgADkooE7pkyyPPY9v8xfM4lgI5KkVLXjMfWFx0gbKjBa985i+e1/LrHn7Na22D9waQv1KEPADpyvoDm/svd7/oa1MvJHLKC5C/UoS8pANLt3dJ3XfrZj1777McGJGreq06Wb4f1y7ZqPM/e37m5WCkPI6FvJ5A/kp9Tl74/JBaB6BOHJy2rn16+Peke0gCecs64OVRDtLH8QnoTqVoyasnKeh1onshOFjOefec2MfSqg8fqEgS94Gdx3skKOOoGRO187ldmZeHdDQS8DjGAJ1ByKuhUwEvSRr1SbrfaBZloRh4LyAOrXdA9yt1d203aeSA2e5BtZn1OYCXU/R1J+AVX9Z/67yHgM6hJa8L5BGogjLq3d0jPIQMV8OVG6892sV7COgsQl6HZBns4rYdhLu5mltz3bmaq0ENDoc9zt0rH17n8qDFHig2Ql7O8g53UiXg1Qt3cevE3dotjDBQPLym5UHAA4qPkJeTNA6ocR/A9bbbSsALWvOC1j/uAgIUSyPHjVo4JgC9gZCXg1YD3pyrT9eBXfta3majAe+K/u1aMTBtxLxgvSDsxR3UOWcP6D1p/MOZdEwI49gAdF7dkPfUU09p1Kj0L8IdM2l86tvsVtPPnt7U8mdcearGHzdR+3fvk8Xs+wdue1iv++x8PXDbw4nbOEfn6KV6qZ7TczXL+tRpz2n5Q5P03Nm1lwv+s1+v9SPmF/V1tD4r7HPrduz7bNU6Ho2eOFoLzx7ZivdSvfSI5R7Vo5qu2se1h1dv0RlXntpaJUto1OhRGnPMOB06cChxmSeffDLHGqEI6oa8HTt2ZBPyjhqb+ja70dAPn9SU06c0tOwJr3/x8O+jxvRp9Pgxsfv+le+ZLUmJ252lWZqu6dqjPTXL++aM3fq9rRO15/Tk5YID+TZtq5SpkWUW9XW0UaMK+9y6Hfs+O/WOR0eNO0o6vfJ78LefdBw5R+foQT2YuK3w8Qz12ehRGj1hjPyQJy6zY8eOHGuEIqgb8s444wz19fWlXvB/fP3O1LfZbZrpFpm/eJ52bz3cmnbUqMnau3NPYnetJP38tp8fMa+Zc/De33DtDoe8WZo1optm6ssmNrGV3jFm0vgRrwfyw77PRq2r7oedI+m2+lfhH6/jNajBI44HYUU9NmSlb8IY7dm2S4f2H0xcZvbs2TnWqLjMbIGkGyT1SVru7ktjlrlA0t9IGiNpu7u/PtdKpoTBkDPU6DkpaZ27EhywBzWYyvbicMJ1sXHFZfHE3ft6YfUraoImNH2hFtBLzKxP0o2SLpE0R9JVZjYnsswUSTdJeqO7v0LS23KvaEq48CJj8xfPq/nB2U7Ai257lVbFnhAdTIcP6o0eyIPASLgrjuA9E/feK/rJ8rWee9HEBbskwfFgszZLqlyMFVgxMC324qxmy48qw2uArnSupM3u/rAkmdltki6XdF9ombdLut3dH5Mkd9+Wey1TQsjLQVLQy+IgF4SxuFAWDntBeKsV9uoFPA7SvaesV0SXqYWylXAXFg54cdNpKVPgxpFO27xJk58+NtVt7t69SZuk2Wa2ITR7mbsvC03PlPR4aHpI0nmRTb1M0hgzWyPpGEk3uPtXU61sTgh5OQkHvbwPanEtftLIsCeNPOCH73qRtE30njK+bkUaxLveP4v1Al7Q2i8dGfCu6N+uayYcSCy7ldY8oAM2ufu1NR63mHnRq11GS3q1pIskTZD0EzNb6+5Hngjf5Qh5OcrzwyVaVjAdF/ai4ub38gcjyq0I792klsikgBcX7gJx4Q4okSFJJ4amT5C0JWaZ7e7+gqQXzOxHks6URMhDvpr9AIsLewHCHdB90gp4cd22BDyU0HpJs8zsFElPSLpSlXPwwr4t6fNmNlrSWFW6c/8611qmhJBXUrXCHsEO6LxGL9iKLhf3z1q0i7aVcEdXLYrA3Q+Y2fWSvqfKECq3uPu9ZnZd9fGb3f1+M7tD0j2SDqkyzMrPOlfr1hHySo5AhzLq9gtQGrlQpNGLSZLOwUtDGlfdd/trgeJx99WSVkfm3RyZ/oykz+RZrywQ8gAURlmuom0n4KXRihdcsNXIPWwbwZW2QDYIeQB6WlLgCQJOWc81TaMFLy7cZTlmJmEPSBchD0DPauhWXRG9ECDqDaJeTy8GvDAGUgbSQcgD0HNaCXe9pt2glyQa3uK6bzsZ8BrBeXxAY7h3LYCeUi/gza1+SfFdtcG9XHvh/L1Wgky7rXjhgDdY/eqmgBfohdcP6DRa8tAT6g0nwX/25ZD0Pgju5BAEm0aDSbe/b9oJMoMabCjohUPd5j07Rqwvcd9qoJcR8tDz+I8ekkYEPFQ0co/qWvIOeI10uYfrFIR0LtgA4hHyAPSM6Id48OFe6zZeUdFhP7q5Na+V8/LC96dtVd5BOdrdHidp2Jbw/ml0AGmgLDgnD12Pljo0ql43bdK4bt38HmslnKTRApdXK17QzR4+lzJO+LFWQmyvnIcJpImQh8LgAF4u0Va8pVpa8z6tUu3g0s0hoNmgF33+gzFfcfZoT26teAurX810JYeDYKutld38OgNpo7sWXS+roSRQDNHg1k5XZbd26TXz/q8XcIPHO3nuYlq3Wmvkjhvh/dEr3fRAWgh5AAqj1l0uGtGNH/ztBLy4/dDKOXvhfdIN/3DN1dy6ITXuOaZ1GzagV9BdC6AnzV88b0T4aDfgBbohxLRqVeSr3nJREzRh+PEk7YbgtFrxkrptg27gYJla5/v18msNNIKQB6AwitZKk2cImau52qM9R8yPC3VptHam1V0ctEzGhbuouZp7RCgk6KHI6K5FYXRbNxuyF/6ADgJevfcBH+rx4kJXFn9T4e7iRgdsrqWR4Vfi1inaPwRAHEIeekLSxRcEO4QV6f2Qdytesy1rSWMW5qmVegfrhINeN56LCaSBkIeeERyEOSAjUOT3QdZXlYdbwIKgtFM7dZfuGi4/K/WGumlGKy2BceGQ4wqKiJCHnsOBGO3olSF5WqljrW7I8D9JUvz5i+fr/BHLNqLdfRnXZXtF//YjlgvfY7cZwbai6wdl0m2LIiPkAUBIUsDJ8/6orQa8JFnVOY2w3EjAa1Yz2+D8PBQZIQ9A6dUKQXm2+rVaVlrDxzQi6/2xYmBabEi7on973da8WuEuaf1GL9hBcZjZAkk3SOqTtNzdlyYsd46ktZJ+092/mWMVU0PIA3CEPFutOqHbnlc7wanZgFevu7rRgY+7pds7jZY/iXPyysLM+iTdKOliSUOS1pvZSne/L2a5T0n6Xv61TA8hD4WR9IHDgbu+bviw7lZZX/BTb9+3G6YaGecueG5Hz5ioHTt3NVx2u3ULgmm0yzZocYsGuLQCHUrtXEmb3f1hSTKz2yRdLum+yHIfkPQvks7Jt3rpIuShEGp90JTpP/RGnmszH8pl2ne1BPss7/3RblmNrl/vPMQ069SMpLCXBd7n+frYQ5P0erV2MU2Sn2mK3iTNNrMNodnL3H1ZaHqmpMdD00OSzgtvx8xmSvoNSW8QIQ+oLfpBkebBtNHAUoawEg4iaW+36PsuSadaOLthf+fZZZ/GoMj1tHp1LnrOJne/tsbjFjPPI9N/I+kj7n7QLG7x3kHIQ+7iQl/aLVBJ63fDh2cW6G5NXx77tJGWslr1SLqiNo33eZp/K90wcHLUoAa56KKchiSdGJo+QdKWyDJnS7qtGvCmSbrUzA64+7fyqWJ6CHmoKwhHjZw71Ow6wbLhn/WE708ZDGga/r3WCeiNlnH+Z88fsWyaJ6Nn0ZJZa/iMuP0RXp7hI46UVwhp56repAsu8g4szVzIESd8m7Mk7XTb0oKHiPWSZpnZKZKekHSlpLeHF3D3U4Lfzewrkr7biwFPIuT1jEY+dNI+uIfLbCasNbNOMyPeJ91wPCyrwNJMt3Ajy7T7WsWVU6+7K+kuA40MoIvDsuzGbCZcxr1mnXi92r3gqdH72KZ1Xh4teOXm7gfM7HpVrprtk3SLu99rZtdVH7+5oxVMGSGvB+TdtZFHeeHWuDylOZ5YrdaweuW0E/TCr89CLRzeh+EPwaTWizRuI1V0nfiHqtFye0Wr+yfcMp8GxsVDHHdfLWl1ZF5suHP3d+VRp6wQ8gqkF8Y2iwt30f/QkwJK0n/yjXbHhAPOIi1qqHu33nZqzQvPT6uFMdw9W2sfhgd+Dd+QvRnd/D7qlCz3SaOnBDSyfp669X0Sd1yI3q8WKDpCHkbIqjUhqeUuaWT7VtQKUnHlh8NPoyGs0RbI6IdJXNhrtjUvLuBlObRE0jmJRdYtYy02O/5cL70+0ecV/ptIq6U5KeDRioeyGdXpCqA7rF2yMZOAt7D6Nbf6Fbiif3vuA5s2ck5fLa0GvOg2wh9kzZ7r10rAq9d6keeFF1m9z9IS/vCfv3je8Hen61JPN+/TqEaeVzstbgQ84DBa8gqo2W7brD8gsrgBeSA4oIc/FJppBUn6MDlqxsQj5jVyFWBUXCBstvUwkBTwBvor43j2D6wbnhf9oGv2nMGwtD4Ue+F0Aql765fn/WmzlPXxJnqPWrpoUWaEvALLY+y5WuqdN9aucMBbpVVHtMLUG14k/N99WLCdpLAYvQow6UMkrZPH1y7ZWDfgJUl6js3c5zQt3RqeulXce6/XA147VgxMa+gYEhfwyrzfUG6EvAJrZwyudkWDVTMBL667JVg/+lhSiAlbpVVapEV110u6x2d4XwWtedFgFx2eJM2AF95eIwEvGn67Wa+07nWDbn8t0xBuLY8bUqVe0GNMPGAkQl5BdWqE+VZCTrTLdVBbj1jmzwYqP5dq+Yj58xfP03zVDgjRwBkNP/UCRvTxVUuSP2yjXbpJrXy1zg9qZIgUaWT3bKDVgJdnyIp7LxL2Rmrl77Uod3SpF/QakWULXhkvSELvIuQVUDceeOJGrM8rkES7VpsJeI2Wn9StK9W/A0dSt3C0u7ueXgt4SVccB3cb6cb3cTfrhf3VaHht5fzXRlrx2tlHvXRxCxAg5HWprWuf1tSXHz083ejFBI0GkDitXAzQrFbPl2nn4By0BqQR8JJEX5+lWlp3mUB0XlJraFxXVXh/RsvM40O/1Y5Ucg4AABHfSURBVDAW/gBPuqCgW4YzQWeET40I/y1E/w7iTuEI1k8D70P0MkJelwk+NGfMf5H27twz4rEsDyrhgNdO2Es6MIdlGbaS6hT+mXWZtc71qzfgbb3u7lrDQ+T9odNsy0bchTBxz3OndrZXsQJodpy8Mmr2/DtCGcqIkNdlgoP7RX//+hHz0thmkuBDN60hGqJBL4+WtDjh553HGFnBtls5Zycp3IXP6av1unSqy7WdcqPPNZiO3p0jz/dML+uFFqdmzxVO6rZNCnh5teIFj3XTvgXiEPK6UJ4HjqzuZRo9yIa7EvN8fp04CDdbZr2A120hp5U6NNMqNUETdLpOHzGvDFeWNiL699rIfunmi1ravZVbpxH00O0IeWhKtFWwXisho8w3Jql7tgj7L+5ii1rjJ14z8UAk4mVTF6l39mu9eyOHJQW/bg8k3dZF3czdaLp5v6LcCHklkHSwaueG9UkDD9daB4eF91fSXTHK3HoVvWtBu5Len936AZ1U30ZuqZfU1V10aXbVNnLMLMI/YCg+Ql4JNPofcq27P7QyD/Hi7kMrFXd0/uj7L8u7oDTbEtSN79u4lk+psbEnw+c19qK4Y1W9v4fofok7t7jZ1zlaj2gLdLB/u/H9A4SN6nQF0B0aDXhI32D1a1X1K043dWO1otn3UiuteEULeK1KCnjd+HwbVSvoBX8/UUEwa+cc0oXVr7nVr6C8TlzNDrSClrwSSzpwcvDKTvQWaeGfZRA817Ru+ybVD0ZpXTXeCc0OClzrauxu7ZqOSup5CD+fuH1Sa9DzVqR9i0KgE2jJK4leOLiXwRG3SGvig6jXW/PiJLXYNdqNmzTWYPC1SIuGW2GCeeF1e2GfNvIeCbdmFaFFuNbxKhrwwq1s0sig22orXrT1LrBUS/XM4ic4nqJnEPJKhANTd5i/eN6I7yRZDW/TKcFzXaVVDZ0zFg16jYyxFv5wjhuHb67mapEWHRH2elWj4S6sV55vvb+PWlptxUvaN4Ma1FIt5RhaEGa2wMweMLPNZrYo5vH/Y2b3VL//08zO7EQ900B3LUbgINZ5SeGuiK9NUotd0MIXjK/YSMBbpCOO1bFlrBiYNhwAgzDQLV2ZtS6SanRgbKn20Crd8DybUWuf1Bo8vNnnmXRxUBDu5qu39hvimVmfpBslXSxpSNJ6M1vp7veFFvuFpNe7+/+Y2SWSlkk6L//ato+Qh2G9dvAvoqK13iWpdZ7TqbtHjzinqpH3ZXS/tXLVbjcMGpwUZhptmYo7j6xXr7RtVSuteNEracOtoxwXC+dcSZvd/WFJMrPbJF0uaTjkuft/hpZfK+mEXGuYInP35AfN/Ac/+EEmBT+z9elMtlsUo8b06dDBQ9Kh5NenXbu27pYkHT1jYmZl9Kq+caN1cO+BEfN2bd2dyb4KXgdJmqzJw7+H7+FalNdo19bderFeLEl6eMLI/Xvqnsr/nLtm7tKuJ3bVfM7RfTZBE2K3KUl7Jh4lSZqw+4UR5UjSHu1JvFdu3vs8/JySHD1jYuJywXsn2BdS5flJjd8P+OgTjtauoV0N1aNZ7fz9xP2NhJ9n4Ek92VL94ra/UztzfQ+MGjtah/YflGp8Jh8340WZlX/hhRfK3S2zAkLM7JNf0pc+9nq9vv7CTfiZfqY36U1fdPdra5T9VkkL3P2a6vRvSTrP3a9PWP7DkmYHy/eaui15r3vd69TX15d6wV/6+I2pb7NIjpo5WXt37tGBXfsyL2vHzvoH9bIZPelo3fXhu2IfS/M/+0bOj5q/eF4hXqPo+IDBXS2CLtmx1f06+cXTdWDnrsTnnDTOYNB6N9B/ZK9K/8A6qRoKgq7geldg1nqd0+jejbYe3aX677dgn0TfN+GWzG/qmy3X6fzPna+7/jC+HuH61Hs/1ntft9qVWuuK10ENNrQPa20/unyef3fHnHqcXnh8RyXoJXjLVW/LrT5Ze0SPaKqmprrNh/SQJM02sw2h2cvcfVloOi7IxiZrM7tQ0nsknZ9aJXNGdy3QpLTO3+qVE+DTMnzhxZJVI7pi0zzXKTngVawYmNbwoNP1Xp92unejoSVJo9tudpzLVt97zYalVrbRyjaD13ShFh6xL5opj67ZQthUqyVPlfPwTgxNnyBpS3QhM3uVpOWSLnH3Z9KtYn4IeSiFNM+3yjPgFfFDJ8/nFA14Uvrj5TUT+pPuZtFsSIu7GKHZMFPv1l2N7qekq5zjttGt7+durRcysV7SLDM7RdITkq6U9PbwAmZ2kqTbJf2Wu/88/yqmh5CHQku7tSytDwPu95uuK/q3j2jFiwt30uGu4STBvm82ADXyT0QzAa8R7bxPkp5f+JzQZtYL1GqZzPJ9PVdzY7vf+VtClLsfMLPrJX1PUp+kW9z9XjO7rvr4zZIWSzpO0k1mJkkH3P3sTtW5HYQ8FFY7Ae+oGRMz/4AIBwo+jFoXnIsXDnbSyHPvpHRapuKmw9tNei3jAl6r5wO2q5HnF8xrZJ8lBbt27h1bz6AGuRMFWubuqyWtjsy7OfT7NZJ68kKLKEIeCinug6xbg1S31qubBRdd/P/27ic0jvMO4/jzWkaRmtKDrRobOYUQB1KBhFlM7IOh7iHFQYH00JqQ0kOJWuI0lB5KKwoRdHvRtYfYxohcE3xoi6lcfNOhuCvZDCqyg1tESrEcmyKrJq4pTU3eHlazenc8s/80M+/s7PdjDNrZ9erd8Y720e/9F7dUSq/hTuosACWFi6Q196LP2a49eb8f3Nc3qlEd0ZGmMW5ScpvbBby0hkdEz2Hc+nhU8YCnEfJQOv0U8NCbuKAVTqqY10Ju3zM8vqa1pspXGPSKPrnGDWmTmtS61htf9yLr/YGje/kmBXmud6COkIfS4wd+OURnpIZVvOjuGD6EQc/VTcDL6z3ayVIkcTpZXDnL7tlQ2O5Wu39wvQM7CHkoNHdAe6cVOn7Il1u4R2slqP9fL71WP+4uxdJtBa2TbdPacffL7aSrOM/3aXRMYNz6gpI0M/pEv4zZBs4VfY1S9lumhdd/WMlLOq9c+0AzQh4Kwx3LFP2QTfrQLcJWVMhHWkuE9PK8cRWvuD1xo4/PuvuyF2Hbetn6Lfo8eW2Z5v5fdrvcDDDICHnwqtMw126ANzNUEdUu6HW6Dl34HmvXtemGpjDwtZq4kNd7Nm5mbzcB70xl86lqnivrKl4rXPNAa4Q8eNHrgPRB3ngd3dvtosG9CoNRXDDsZMmVLHVbwWsV8Hwh3AGdIeQhd90EvG4HiANxegkFrSqBnQQltwIWN44tL+Hr2BmL13nIiwa8uGplHgh1QG8Ieei5Syvt7xXVbjPyPGbzof8N+nsjes21Cmfr/3nYMtglGfRzDBQVIW+AxW2S7lYcFrXY9AGR5w/ybgIekKfdTFiITsjI85oKZyWHOqmSt3utXIdAsRHyBlAna2WFy0G4S0GkMX6o2wVi3e2Lkrq6qCKgqNpNWshD0vIwYUCLLjDciXDhaQDFRsgriGjwmfjBhO6v3NfWra3UQkz4PWY123S81W/rl4J6wGoEvuruqxDt/l24HpZUD6Ksag9f3PXZpPp1EIY297oJKscTn+OM6nvqXgrGGgs2+3zvxoW+8PWdtCcbv1i12jKuCK8DQHt7fDcAda1+WO5ma6RadbXxV+ou4EXvj1b8stqyyT0Xi86fpMcAWQrfa4tajK1eNQW8yvnmv9vCgHdi7mjh3rtJbYpWIHe2jSPgAf2CSl6BJC0EHB7r9IdqUvhyA16n44riVrsPq3pZjsfph30/gUbAcwLdzp1nVQmW+6Zr0w2z4S904fXv7glMuAP6ByGvgE7MHdVndz5/6njcel+tdolwdTvmJhTdHzQ0qcnUxukl6eR1AXlyd4toF/CknerXohb7Ihw9e+hLkpr3AQ7b7W4bB6A/EPIK6tCJr2rf17+sP/18Jfb+aMDLm1uZyLod0apeP3xYYoBsBzpXJaiPw5vXQr07tI8CEtcXUB6EvIKLVrPc2+4P43bVrujWTJeCsZ67bEPuRuF5Bz3AB3e2dxjkourVu3t9U70DUF6EvD4R/bDo5cMjLujFSeqi9SkacoE8RWfZhrPOXWG37BLdmwAKgpBXEq3Gr7n3hTNVW62R12m4y2vHCfc1+eyixmBrTEyoJs/2JtgBKBJCXsm0CkDufQ90t6kyIT09OSMaAPNar65dtyzj8+AT7zmgvxljTkv6jaQhSQvW2vnI/d+T9Ivtm/+WdNZa+5d8W5kOQt4Ai1b/ouEtrPglLZVShA87KnsAgE4ZY4YkvS/pFUkbkq4bYy5baz92HvZ3Sd+w1v7LGPOqpIuSklc8LzBCHhK3PZLi96bMq2sWAICUvSxp3Vr7iSQZYz6S9LqkRsiz1l5zHl+TdDjXFqaIkIenJM3gzRoBDwCQsXFJd5zbG2pdpXtL0h8zbVGGCHlI1A8Bj65aAOg/K1rRPd1L9Tk3tSlJLxljbjiHL1prLzq3Tcw/tXHPZ4z5puoh72RqjcwZIQ8AAJTFbWvtj1rcvyHpOef2YUmfRh9kjJmStCDpVWvtg3SbmJ89vhsASL1V5KjiAeg3teoqQ1P8ui7pRWPM88aYYUlvSLrsPsAY8zVJv5X0fWvt3zy0MTVU8lAYne5VS7gD0G9q1VVNa1rTGtekJrVWXdODubu+mzVwrLVPjDHvSrqq+hIqH1hrbxlj3t6+/4KkOUn7JZ0zxkjSE2vtMV9t3g1CHgqHEAegTE5VT+tc5ZiCxpFNXQomNV9l6zsfrLVXJF2JHLvgfD0jaSbvdmWB7loAADKwvzquU9XTja0ig0p9Eud7wT2taa0R8PZXx721EeVGyAMAICNhwJOkSrAsSfp15VCjq/bcH44lLjgP7BbdtQAAZKC+NeROyAsreTMLP1QtkILKjCrBgqfWYRAQ8gAAyEBTN23lvBQcVyVY1qxmdWTmBUlSJVhgXB4yQ8gDACBl+6vjCma2w50CKTjbuO/IzAtS5bzW3/mWlgh4yBAhDwCAFNWqq5rV6fqN7XAXjscLu2wrCxXNa14nRMhDdgh5AACkaFrTOlPZVCXYbDreCHjBMt20yAUhDwCAFD2Yu6t3dFf7q/WFj89UNhsB73cLv6/PrH2NgIfsEfIAAEhZrbqqoHJM4ezasHr3gOodcsQ6eQAAZOhSMEb3LLwg5AEAkLJpTUuqB7x5zRPw4AXdtQAApCgci/desKaluRvMoIU3hDwAAFJyqlpfOmVNa42tywBfCHkAAKSgVl1VTauSRPcsCoGQBwBACgh2KBomXgAAAJSQt0reW7/6sa9v3Rdu3rypgwcPamxszHdTBtLS0pK+8+Z3fTdjIHHu/eHc+7O8vKypM1MaHR313RSUCJU8AACAEiLkAQAAlBAhDwAADAxjzGljzF+NMevGmNmY+18yxvzZGPNfY8zPfLQxLcyuBQAAA8EYMyTpfUmvSNqQdN0Yc9la+7HzsC1JP5H0bQ9NTBWVPAAAMChelrRurf3EWvu5pI8kve4+wFr7T2vtdUn/89HANFHJAwAAebp5VVdXJX2R8vMaSc8aY244xy5aay86t8cl3XFub0g6nnI7CoOQBwAAcmOt/VDSh56+vYk5ZnNvRU7orgUAAINiQ9Jzzu3Dkj711JbMta3kraysyJi44IssPXr0SFtbWxoeHvbdlIH0+PFj1Wo1380YSJx7fzj3/jx8+FBBEGhoaMh3U8ruuqQXjTHPS7or6Q1Jb/ptUnbahrypqSnedB7cvn1bBw4c0L59+3w3ZSBdu3ZNR4+yD6UPnHt/OPf+BEGgiYkJjYyM+G5KqVlrnxhj3pV0VdKQpA+stbeMMW9v33/BGHNQ0g1JX5H0hTHmp5ImrLWfeWt4j9qGvJGREUKeB3v37tXw8DAXvCd79uzh3HvCufeHc+/P0NCQnnnmGc5/Dqy1VyRdiRy74Hx9X/Vu3L7HmDwAAIASIuQBAACUECEPAACghAh5AAAAJUTIAwAAKCFCHgAAQAkR8gAAAEqIkAcAAFBChDwAAIASIuQBAACUECEPAACghNruXbt3b9uHAACA/vYP3w1A+oy11ncbAAAAkDK6awEAAEqIkAcAAFBChDwAAIASIuQBAACUECEPAACghP4PGRXZjWGe4swAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 720x720 with 2 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig = pl.figure(figsize=(10, 10))\n",
"ax = fig.add_subplot(111, projection=albers)\n",
"\n",
"cbar_kwargs={'ticks': clevs, \n",
" 'norm': norm, \n",
" 'extend': 'max',\n",
" 'fraction': 0.040}\n",
"\n",
"bom.precipitation.plot(ax=ax, transform=aeqd, \n",
" cmap=cmap, \n",
" norm=norm, \n",
" add_colorbar=True, \n",
" cbar_kwargs=cbar_kwargs)\n",
"\n",
"ax.gridlines(crs=ccrs.PlateCarree())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Easily read multiple consecutive datasets into one xarray Dataset (needs dask)\n",
"\n",
"This simplifies the workflow at least for netcdf source data. But this can be done for other source too, but will take some tweaking."
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [],
"source": [
"filename = '/home/data/kai/projects/pysteps-data/radar/bom/prcp-cscn/2/2018/06/16/2_20180616_10*.prcp-cscn.nc'"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<xarray.Dataset>\n",
"Dimensions: (time: 10, x: 512, y: 512)\n",
"Coordinates:\n",
" * x (x) float32 -128.0 -127.5 -127.0 -126.5 ... 126.5 127.0 127.5\n",
" * y (y) float32 128.0 127.5 127.0 126.5 ... -126.5 -127.0 -127.5\n",
"Dimensions without coordinates: time\n",
"Data variables:\n",
" valid_time (time) datetime64[ns] 2018-06-16T10:00:00 ... 2018-06-16T10:54:00\n",
" start_time (time) datetime64[ns] 2018-06-16T09:54:00 ... 2018-06-16T10:48:00\n",
" proj (time) int8 -127 -127 -127 -127 -127 -127 -127 -127 -127 -127\n",
" precipitation (time, y, x) float32 dask.array<shape=(10, 512, 512), chunksize=(1, 512, 512)>\n",
"Attributes:\n",
" Conventions: CF-1.6\n",
" title: Bias Corrected Radar Accumulation\n",
" institution: Commonwealth of Australia, Bureau of Meteorology (ABN 92 6...\n",
" licence: http://www.bom.gov.au/other/copyright.shtml\n",
" source: rainfields 3.0.28 ho-rainfields 2018-03-23\n",
" station_id: 2\n",
" station_name: Melb"
]
},
"execution_count": 14,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bom = xr.open_mfdataset(filename, concat_dim='time')\n",
"bom"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<xarray.DataArray 'precipitation' (time: 10, y: 512, x: 512)>\n",
"dask.array<shape=(10, 512, 512), dtype=float32, chunksize=(1, 512, 512)>\n",
"Coordinates:\n",
" * x (x) float32 -128.0 -127.5 -127.0 -126.5 ... 126.0 126.5 127.0 127.5\n",
" * y (y) float32 128.0 127.5 127.0 126.5 ... -126.0 -126.5 -127.0 -127.5\n",
"Dimensions without coordinates: time\n",
"Attributes:\n",
" grid_mapping: proj\n",
" long_name: Accumulated precipitation\n",
" standard_name: precipitation_amount\n",
" units: kg m-2"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"bom.precipitation"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Display the source blob
Display the rendered blob
Raw
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
Display the source blob
Display the rendered blob
Raw
Loading
Sorry, something went wrong. Reload?
Sorry, we cannot display this file.
Sorry, this file is invalid so it cannot be displayed.
@cvelascof
Copy link

in xarray_pysteps.ipynb [14], by using

bom = xr.open_mfdataset(filename, concat_dim='valid_time')

all files are correctly concatenated by the end time of the accumulation period.

@kmuehlbauer
Copy link
Author

@cvelascof Yes, you're right, this would work too, thanks! Regarding accumulation period, how is it with the MCH GIF files? Is it the minutes before or after the timestamp in the filename?

@cvelascof
Copy link

Not sure on MCH GIF files. Best assumption is always minutes until timestamp. Perhaps @dnerini knows.

@dnerini
Copy link

dnerini commented Jun 26, 2019

It is the accumulation in minutes before timestamp for all the MCH QPE products.

ps: I really like the way your're putting together xarrays and pystpes, very impressive!

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment