{
"cells": [
{
"cell_type": "markdown",
"id": "360075b1-2ac9-4ca1-9615-403e5481d250",
"metadata": {},
"source": [
"# Calculating zonal means"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "ec783093-5325-45d2-8712-a741e8faa220",
"metadata": {},
"outputs": [],
"source": [
"import intake\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"\n",
"from gridlocator import merge_grid"
]
},
{
"cell_type": "markdown",
"id": "de341c08-79f1-4a4c-a93e-997fbf8e87cc",
"metadata": {},
"source": [
"In this example we will compute the zonal-mean air temperature.\n",
"We retrieve the 2m-air-temperature from the `dpp0067` NextGEMS simulation using `intake-esm`."
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "06aee1d5-5f60-4586-9b10-87a5315164b3",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"--> The keys in the returned dictionary of datasets are constructed as follows:\n",
"\t'project.institution_id.source_id.experiment_id.simulation_id.realm.frequency.time_reduction.grid_label.level_type'\n"
]
},
{
"data": {
"text/html": [
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/html": [
"\n",
"
\n",
" \n",
" 100.00% [1/1 00:00<00:00]\n",
"
\n",
" "
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"catalog_file = \"/work/ka1081/Catalogs/dyamond-nextgems.json\"\n",
"col = intake.open_esm_datastore(catalog_file)\n",
"cat = col.search(\n",
" variable_id=\"tas\",\n",
" project=\"NextGEMS\",\n",
" simulation_id=\"dpp0067\",\n",
")\n",
"cat_dict = cat.to_dataset_dict(cdf_kwargs={\"chunks\": {\"time\": 1}})\n",
"\n",
"ds = merge_grid(list(cat_dict.values())[0]) # Include the grid information!"
]
},
{
"cell_type": "markdown",
"id": "d844c228-de01-4188-8a39-e082385a8843",
"metadata": {},
"source": [
"## The general idea\n",
"behind the zonal average is to calculate a weighted average of values in a certain latitude bin.\n",
"Therefore, in a first step, we count how many cells are in given equidistant latitude bins."
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "a8430e1a-0f1b-4148-b15f-4dac307fbe86",
"metadata": {},
"outputs": [],
"source": [
"hist_opts = dict(bins=128, range=(-np.pi / 2, np.pi / 2))\n",
"cells_per_bin, lat_bins = np.histogram(ds.clat, **hist_opts)"
]
},
{
"cell_type": "markdown",
"id": "a10c6bda-54c4-491c-9925-5b308ddfb3ff",
"metadata": {},
"source": [
"Now comes the trick! In a next step, we will repeat the histogram but account a weight to each cell.\n",
"Usually, `histogram` will weight each data point with one, i.e. it will count the values in a certain bin.\n",
"Here, we will weight each data point with the 2m-temperature. Thereby, we will compute the cumulative sum of all temperatures in a given latitude bin.\n",
"\n",
".. tip::\n",
" The `np.histogram` function is more efficient when passing a range and a number of bins.\n",
" This is because, when constructing the bins internally, the function can assume equidistant bin sizes.\n",
" This is not the case when passing a sequence of bins directly."
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "5bbb275b-3389-4008-b37b-22fb37b942c1",
"metadata": {},
"outputs": [],
"source": [
"varsum_per_bin, _ = np.histogram(\n",
" ds.clat, weights=ds.tas.isel(time=1, height_2=0), **hist_opts\n",
")"
]
},
{
"cell_type": "markdown",
"id": "5d43cca2-88b3-4d44-8d7c-08d65ffd0a07",
"metadata": {},
"source": [
"The zonal mean can now be computed by dividing the cumulative values of the temperature with the number of cells in each bin."
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "70c53f39-1422-48eb-9a78-aad2b8f83d9f",
"metadata": {},
"outputs": [],
"source": [
"zonal_mean = varsum_per_bin / cells_per_bin"
]
},
{
"cell_type": "markdown",
"id": "9ad77d8e-2bea-4bb2-a813-c21e2efadae3",
"metadata": {},
"source": [
"We can check our result by plotting the zonal mean as a function of the latitdue bins.\n",
"While doing so, we will scale the bins by their area so that the visual appearance of each latitude bin represents their actual proportion."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "718e97e5-2ada-480b-b4d7-3df5e41e705e",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"Text(0.5, 0, 'latitude')"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEKCAYAAAD9xUlFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvtElEQVR4nO3deXxV9Z3/8dcnO9kIhARCFhL2fY2gUFCsK+4LiFrr0inVsa3btGrtdGw7dmw7dX5Oq2Np1Vo3REHFui8oymoSdsISIISQQEIgGyH75/fHPaQRAgTIzcm99/N8PPLg5pxzT943Xu8n3+/5nu9XVBVjjDEGIMjtAMYYY7oOKwrGGGNaWFEwxhjTwoqCMcaYFlYUjDHGtLCiYIwxpoXXioKIRIjIKhFZKyIbReSXzvaeIvKxiGxz/u3hbE8XkcMissb5esZb2YwxxrRNvHWfgogIEKWq1SISCnwF3ANcCxxQ1cdF5CGgh6o+KCLpwD9UdaRXAhljjDkpr7UU1KPa+TbU+VLgKuAFZ/sLwNXeymCMMebUhHjz5CISDGQDA4GnVHWliPRW1WIAVS0WkcRWT8kQkdVAJfBzVf2yjXPOAeYAREVFTRg6dKg3X4Ixxvid7Ozs/aqa0NY+r3UffeOHiMQBbwI/Ar5S1bhW+w6qag8RCQeiVbVMRCYAbwEjVLXyeOfNzMzUrKwsr2Y3xhh/IyLZqprZ1r5OGX2kquXA58AlwD4RSXKCJQElzjF1qlrmPM4GtgODOyOfMcYYD2+OPkpwWgiISDfgAmAzsAi41TnsVuDtVscHO4/7A4OAHd7KZ4wx5ljevKaQBLzgfNAHAfNV9R8ishyYLyLfAwqAmc7x04BfiUgj0ATcqaoHvJjPGGPMUbxWFFR1HTCuje1lwLfb2L4AWOCtPMYYY07O7mg2xhjTwoqCMcaYFlYUjDHGtLCiYIwxpoUVBWOMMS2sKBhjjGlhRcEYY0wLKwrGGGNaWFEwxhjTwoqCMcaYFlYUjDHGtLCiYIwxpoUVBWOMMS2sKBhjjGlhRcEYY0wLKwrGGGNaWFEwxhjTwoqCMcaYFlYUjDHGtLCiYIwxpoUVBWOMMS2sKBhjjGkR4q0Ti0gEsAQId37OG6r6HyLSE3gNSAfygVmqetB5zsPA94Am4Meq+qG38pnAVV3XSEFZDQUHaqg83EBdYxN1jc2er4Z/Pq5taCIiNJgRfWMZmdydgYnRhAbb31HGv3mtKAB1wPmqWi0iocBXIvI+cC3wqao+LiIPAQ8BD4rIcGA2MALoC3wiIoNVtcmLGY2fq65rJHvXQVbuKCMr/yDbS6spO1R/3ONFIDwkiPCQYCJCg6iqbaSm3vMWDAsJYkBCNL2iw+gZFUaPyDASYsIZ2ieG0SlxJMSEH/e8VbUN5JVUs62kmrLqehqbmmloVhqbmokKD2F431hG9u1+wnMY0xm8VhRUVYFq59tQ50uBq4DznO0vAJ8DDzrb56lqHbBTRPKAicByb2U0/qXicANb9laRW1zJ5r2VbCzyfDU1K8FBwqjk7lw4vDdp8ZH06xlFWs9IekSFthSA8JBgQoMFEWk5Z1Ozkl92iA17KthYVEleSTUHDtVTcKCGA4fqqaptbDk2qXsEI/rGEhYSRP2RlkdjM4UHaiiqqG0zc0iQ0NisLd/3jg1nVHIcZ/fvyZSBvRjSO4agIGnzucZ4gzdbCohIMJANDASeUtWVItJbVYsBVLVYRBKdw5OBFa2eXuhsO/qcc4A5AGlpad6Mb3xAxeEGPtywl7fX7mHZ9jLU+XyNiwxlWJ9Y7jp3AJP692R8Wg+iwk/97R4cJAxIiGZAQjRXjT3m7cihukY2FlWyrrCc9XsqyC32FKGwkGDCQ4IICwliYkZPBvWOYVBiNIN7x5AYG05ocBAhQZ4CVFnbwKaiypbCs7rgIJ/k7gMgPiqMyQN7MfusVCYPiP9GwTLGG7xaFJyun7EiEge8KSIjT3B4W+92PWaD6lxgLkBmZuYx+43/O3ionsVbSvhw414WbymlvrGZ9PhIfjh9IOP79WBYn1h6x4Z3ygdoVHgIEzN6MjGj52mfIzYilLP7x3N2//iWbXvKD7Msbz/Lt5fx+dZS3llbxODe0dw2OYNrxiXTLSy4I+IbcwyvFoUjVLVcRD4HLgH2iUiS00pIAkqcwwqB1FZPSwGKOiOf6ZqampWy6jr2Vtayt6KW7aWH+GzzPrJ3HaRZITEmnJsnpXH12GRGp3T3q7+ik+O6MTMzlZmZqdQ2NPHO2iKeX5rPz95cz28/2Mz3p2Zw+5SM02r9GHMiouqdP7ZFJAFocApCN+Aj4LfAuUBZqwvNPVX1pyIyAngFz3WEvsCnwKATXWjOzMzUrKwsr+Q3nauhqZn1eyrYsKeC9YUVrN9TQV5J9Tf62wGGJ8VywbBEvj2sN6OSuwdUf7uq8nX+QeYu2c4nuSX0ig7j7ukDuWlSGuEh1nIw7Sci2aqa2eY+LxaF0XguJAfjuR9ivqr+SkTigflAGlAAzFTVA85zHgHuABqBe1X1/RP9DCsKvq3icAOfbynh09wSPt9SQqVz0TY+KoyRyd0ZlhRLclwEvWMj6NM9guS4bsRH2+gcgJyCg/z+gy0s31FGclw3vj81g+smpBATEep2NOMDXCkKncGKgu9QVbbuq2bt7nLWFJazdnc5m/dW0dSsxEeFMX1oIucPTWRsahxJ3SP8qivIW1SVpXllPPHxFnIKyokKC+b6CSncck46AxOj3Y5nujArCsYVBw/V82Xefr7YUsqSbaWUVtUBEBMRwtjUOMamxnHeEE8hCA6gbiBvWLu7nBeW5fOPdcXUNzXz7aGJ3HPBIEanxLkdzXRBVhRMp9pVdoinF29nQU4hjc1KXGQoUwclMHVQLzL79SA9PiqgrgV0pv3Vdby8ooDnlu6k4nCDFQfTJisKplNsL63mqc/yeHttEcFBwuyzUrlmXDKjU6wl0Nmqahv4+/Jd/OXLHZTXNHDZqCQevXKE3TFtACsKxsuampU/L9nO/3y8lZCgIG6elMacaf1JjI1wO1rAq6pt4NmvdvL04u1EhQfz6JUjuHJMX7tmE+CsKBiv2X2ghgfmr2VV/gEuG5XEL68aQS8bIdTl5JVU8W+vr2PN7nIuHN6bx64eaUU7gJ2oKNiUj+a0LVpbxKVPfsmm4kqemDWGP900zgpCFzUwMYYFd03mkRnDWLK1lOufWc7BE0wMaAKXFQVzWl77uoAfv7qaoX1ieP+eqVw7PsW6JLq44CDh+9P688r3z2ZvRS0/fDWHxqZmt2OZLsaKgjll87/ezUML1zNtcAIv/cskUntGuh3JnIIJ/Xrw2DUjWZpXxn++m+t2HNPF2MQp5pS8nrWbBxeuY+qgBObeMoGIUJtewRfNzEwlt7iK55buZHhSLLPOSj35k0xAsJaCabfXvi7gpwvW8a2Bvawg+IGfzRjK1EG9+PlbG8jedcDtOKaLsKJgTkpVeWpxHg8uWM+0QQn85buZVhD8QEhwEH+8cRxJcRF874Us1hdWuB3JdAFWFMwJNTcrv3xnE7//cAvXjEu2guBn4iLDePGOSUSHh3DTX1eQveug25GMy6womOOqa2ziR/NW87dl+Xx/agZ/mDmGsBB7y/ibtPhI5v/gHOKjwvjusytZuaPM7UjGRfZ/uGlTXWMTd72Uw7vrivnZjKE8ctlwm6/Ij/WN68b8H5xDUlw3bn1+FV9t2+92JOMSKwrmGLUNTfzgxWw+21zCY9eMZM60AW5HMp0gMTaCeXPOJj0+ijkvZrGpqNLtSMYFVhTMN9Q2NDHnxWw+31LKf107ipsn9XM7kulEvaLD+fsdE+neLZR/eeHrlunOTeCwomBa1Dc2M+fFbL7cVsrvrhvNjRPT3I5kXJAYG8FfvpvJgZp6fvBiFrUNx10R1/ghKwoG8Aw7/cXbG1iytZTHrx1lNzMFuJHJ3fmfWWPJKSjnZwvX48sTZ5pTY0XBAPDsVzuZ9/Vu7p4+gBvOshaCgUtHJXH/hYNZuHoPf16yw+04ppNYUTB8mruPx97L5dKRfXjgwiFuxzFdyI/OH8hlo5L4/YdbWF1g9zAEAisKAW7z3kp+/OpqRvSN5Q+zxtiwU/MNIsJ/XTeKPrER3PfaGg7VNbodyXiZFYUAVtvQxN0v5xAVHsJfv3sWkWE2P6I5VmxEKE/MGsOuAzX8+h+b3I5jvMxrRUFEUkVksYjkishGEbnH2T5GRJaLyHoReUdEYp3t6SJyWETWOF/PeCub8fjtB5vZXnqIJ2aNpU93W4XLHN+k/vHcee4A5n29mw837nU7jvEib7YUGoEHVHUYcDZwt4gMB/4KPKSqo4A3gZ+0es52VR3rfN3pxWwBb2nefp5fms9tk9P51qBebscxPuC+CwYzom8sDy1YR0llrdtxjJd4rSioarGq5jiPq4BcIBkYAixxDvsYuM5bGUzbKg438G+vr6V/QhQPXjLU7TjGR4SFBPHk7LHU1DfxsA1T9Vudck1BRNKBccBKYANwpbNrJtB6QHyGiKwWkS9EZOpxzjVHRLJEJKu0tNSbsf3WLxdtpKSqjv+ZNZZuYTbjqWm/gYkx/OTiIXy6uYSPNu1zO47xAq8XBRGJBhYA96pqJXAHnq6kbCAGOLJ6eDGQpqrjgPuBV45cb2hNVeeqaqaqZiYkJHg7vt95e80eFq7ew93TBzImNc7tOMYH3TY5naF9Yvjloo3U1NtoJH/j1aIgIqF4CsLLqroQQFU3q+pFqjoBeBXY7myvU9Uy53G2s32wN/MFmrySKh5euJ6z0nvwo/MHuh3H+KiQ4CB+ffVIiipq+eNneW7HMR3Mm6OPBHgWyFXVJ1ptT3T+DQJ+DjzjfJ8gIsHO4/7AIMBuo+wgh+oaufOlHCLDgvnTTeMJDbbRyOb0nZXek+snpPCXJTvIK6lyO47pQN78ZJgC3AKc32qY6QzgRhHZCmwGioDnneOnAetEZC3wBnCnqtrCsR1AVXnkzfVsL63mydnj6B1rw0/NmXv40qFEhYfw729ttIvOfsRrdyup6lfA8W6PfbKN4xfg6WoyHeyVVQW8taaI+y8czJSBNvzUdIz46HB+cvEQfv7WBhatLeKqscluRzIdwPoQ/Fz2roM8umgj5w5O4IfT7TqC6Vg3TkxjTEp3/vPdXCprG9yOYzqAFQU/tq+ylrteyiapezeenD3W5jUyHS44SPj11SPZX13HEx9tdTuO6QBWFPyUZ43lbKrrGpn73QnERYa5Hcn4qdEpcXxnUj/+vjyfjUUVbscxZ8iKgp96dNEmcgrK+f31Yxja55jbPYzpUP920RB6RIbx87c20NxsF519mRUFP/T2mj28uqqAu84bwGWjk9yOYwJA98hQfjZjGKsLynk9e7fbccwZsKLgZ4orDvPvb21gfFocD1xo9/6ZznPt+GQmpvfk8fc3c/BQ/cmfYLokKwp+RFX56RvraGhSnpg1lhC7Qc10IhHhV1ePoLK2kd99uNntOOY02aeGH3lpxS6+3LafRy4bRnqvKLfjmAA0tE8sd0xJ59VVu8mx5Tt9khUFP7GjtJrH3svl3MEJ3Dwpze04JoDdc8FgeseG8/M3N9DY1Ox2HHOKrCj4geZm5SdvrCM8JJjfXT8az7RTxrgjOjyEX1w+gk3Flby0YpfbccwpsqLgB15csYvsXQf5xeXDbV4j0yXMGNWHqYN68YePttoqbT7GioKPKzxYw+8+2My0wQlcO97mnjFdg4jwyytHUNfYzG/ey3U7jjkFVhR8mGf20w0o8JtrRlq3kelS+idE84Nz+/PWmiKW5u13O45pJysKPuytNXv4YmspP714CCk9It2OY8wx7p4+kP69orjvtTXsr65zO45pBysKPqqipoFfvbOJ8Wlx3HJOuttxjGlTRKhnUafyww3c99oamwLDB1hR8FF/WryN8sMN/PrqkQTb7KemCxveN5ZHrxjBl9v28/TntnxnV2dFwQftKjvE35blM3NCCiP6dnc7jjEndePEVK4c05cnPt7Kyh1lbscxJ2BFwQc9/v5mQoODeOCiIW5HMaZdRITfXDuKfvFR/OjV1ZTZ9YUuy4qCj1m18wDvb9jLnecOsHsSjE+JDg/hqZvGU17TwIML1tu6zl2UFQUf0tys/Oe7m+gTG8H3p/Z3O44xp2x431h+eskQPsndx7yvbYrtrsiKgg/5OHcf6wor+MnFQ+gWFux2HGNOyx1TMvjWwF786p1N7Nx/yO045ihWFHzI80t3khzXjavG9nU7ijGnLShI+O+ZYwgLCeLeeatpsEnzuhSvFQURSRWRxSKSKyIbReQeZ/sYEVkuIutF5B0RiW31nIdFJE9EtojIxd7K5otyiytZseMAt5zTz9ZJMD6vT/cI/uvaUawtrOCPn25zO45pxZufLo3AA6o6DDgbuFtEhgN/BR5S1VHAm8BPAJx9s4ERwCXA0yJifSSOvy3NJyI0iNlnpbodxZgOMWNUEteNT+Gpz7dTeLDG7TjG4bWioKrFqprjPK4CcoFkYAiwxDnsY+A65/FVwDxVrVPVnUAeMNFb+XzJgUP1vLVmD9eMSyEuMsztOMZ0mPsv8iwZ+6JNsd1ldEo/hIikA+OAlcAG4Epn10zgyJ++yUDr4QiFzraAN+/rAuoam7ltcrrbUYzpUMlx3bhoeG/mrdrN4fomt+MYOqEoiEg0sAC4V1UrgTvwdCVlAzHAkRW+25qr4ZiBzCIyR0SyRCSrtLTUW7G7jIamZl5cvovJA+IZ0ifG7TjGdLjbJqdTcbiBt9fscTuK4QRFQUSOO3+CiJzVnpOLSCiegvCyqi4EUNXNqnqRqk4AXgW2O4cX8s9WA0AKUHT0OVV1rqpmqmpmQkJCe2L4tI827qO4opbbp2S4HcUYr5iY0ZOhfWL427J8u6GtCzhRS+FTEelx9EYRuQhYeLITi2dy/2eBXFV9otX2ROffIODnwDPOrkXAbBEJF5EMYBCwqr0vxF/9bdlOUnt24/yhiW5HMcYrRITbp6SzeW8VK3cecDtOwDtRUfgzsFhEWv4cF5GbnO2XtePcU4BbgPNFZI3zNQO4UUS2ApvxtASeB1DVjcB8YBPwAXC3qgZ0J+OGPRV8nX+QW89Jt5lQjV+7amwycZGh/G1pvttRAl7I8Xao6l9EpBb4zGkd3ADcCUxX1fyTnVhVv6Lt6wQATx7nOY8Bj53s3IHi+aX5RIYFMzPThqEa/xYRGszss9KYu2Q7e8oPkxzXze1IAeuEF5pV9UXgV8Bq4CZgSnsKgjlz+6vreGdtEdeNT6F7t1C34xjjdd85Ow2AF5fb8FQ3nehC83oRWQf8AogE4vF0Jx3Zbrzo1ZUF1Dc1c+vkfm5HMaZTpPSI5PyhvVmYU2grtLnouN1HwOWdlsJ8Q31jMy+u2MXUQb0YmGjDUE3guGJMEp/k7iNr10EmZvR0O05AOtE1BWvDueT9DcWUVNXx2+tGux3FmE717WG9CQsJ4r31xVYUXGIzq3VBL68ooF98JOcO9v/7MIxpLTo8hPMGJ/D+hmLrQnKJFYUuZkdpNavyD3DDWakE2TBUE4AuG53Evso6sgsOuh0lIJ3oQvNcEblGRKxTuxPNzyokOEi4fnyK21GMccWRLqR31xW7HSUgnail8BwwBnhPRD4VkQdFZEwn5QpIDU3NLMgpZPqQBBJt/WUToKLDQzjXupBcc9yioKorVPVRVZ0KzAIKgAdEZLWIPCciszotZYBYvLmE0qo6bjgrze0oxrjqslGeLqQc60LqdO26pqCqZar6qqp+V1XHAU/hmZvIdKD5WbtJiAln+hC7wGwC27eHJXq6kNZbF1JnO60Lzaqa7UxJYTpISWUti7eUct34FFtu0wS8mIhQpg1K4P31e60LqZPZp08X8UZOIU3NyqxMu8BsDMDlo5PYW1nL6t3WhdSZrCh0AarK61mFTEzvSf+EaLfjGNMlfHtYIhGhQSzIscV3OtNJi4KIzDwyLFVEfi4iC0VkvPejBY5VOw+wc/8hZp1ls6Eac0RMRCgzRiWxaE0Rh+oa3Y4TMNrTUvh3Va0SkW8BFwMvAP/n3ViB5bWs3cSEhzBjVB+3oxjTpdw0MY3quka7Z6ETtacoHFno5jLg/1T1bSDMe5ECS2VtA++tL+aKsX2JDDvR/ITGBJ4J/XowKDGaV1YVuB0lYLSnKOwRkT/juVfhPREJb+fzTDu8s7aI2oZmbrCFdIw5hogwe2Iaa3aXk1tc6XacgNCeD/dZwIfAJapaDvQEfuLNUIFk/te7GdonhtEp3d2OYkyXdO24ZMJCgphnrYVOcdKioKo1qroQqBCRNCAUz/rK5gxt3lvJ2sIKZmWmImKT3xnTlh5RYcwY2YeFq/dwuD6gl23vFO0ZfXSliGwDdgJfOP++7+1ggeC1r3cTFhzENeOS3Y5iTJd248Q0qmob7Q7nTtCe7qNfA2cDW1U1A7gAWOrVVAGgoamZRWuKuHB4b3pE2XV7Y05kYkZP+idEWRdSJ2hPUWhQ1TIgSESCVHUxMNa7sfzfl9tKKTtUb60EY9pBRLhpYhpZuw6SlX/A7Th+rT1FoVxEooElwMsi8iRw0jtJRCRVRBaLSK6IbBSRe5ztY0VkhYisEZEsEZnobE8XkcPO9jUi8syZvLCu7s3VRfSIDGWara5mTLvcPKkfvWPD+c17uajafEje0p6icBVQA9wHfABsBy5vx/MagQdUdRie7qe7RWQ48Dvgl6o6FviF8/0R21V1rPN1Z/tfhm+pqm3go417uXx0X8JCbHSvMe3RLSyY+y8cTE5BOR9s2Ot2HL/Vnk+kX6hqs6o2quoLqvq/wIMne5KqFqtqjvO4CsgFkgEFYp3DugNFpxfdd324cR91jc1cM966jow5FdeNT2Fw72h++8FmGpqa3Y7jl9pTFC5sY9ulp/JDRCQdGAesBO4Ffi8iu4H/Bh5udWiGs4jPFyIy9VR+hi95c3Uh/eIjGZca53YUY3xKSHAQD106lPyyGl61i85ecaI1mu8SkfXAEBFZ1+prJ7CuvT/AuR6xALhXVSuBu4D7VDUVT5fUs86hxUCas4jP/cArIhLbxvnmONciskpLS9sbo8vYW1HLsu1lXD022e5NMOY0TB+SyNn9e/LkJ9uoqm1wO47fOVFL4RXgCmCR8++Rrwmq+p32nFxEQvEUhJedG+AAbgWOPH4dmAigqnXOKCdUNRvPtYvBR59TVeeqaqaqZiYk+N5F2kVr96AKV9uoI2NOi4jwsxnDKDtUz5+/2OF2HL9zojWaK1Q1X1VvVNVdrb7aNR5MPH8GPwvkquoTrXYVAec6j88HtjnHJ4hIsPO4P57lPv3qv3hzszJv1W7GpcWR0SvK7TjG+KzRKXFcPjqJ55buZH91ndtx/Io3h75MAW4Bzm81zHQG8H3gDyKyFvgNMMc5fhqwztn+BnBnewuQr/hiayk79h/itsnpbkcxxufdd+FgahuaeHrxdrej+BWvzdWsql8Bx+s0n9DG8QvwdDX5reeW7qR3bDiXjkxyO4oxPm9AQjTXT0jhpRW7+JepGfSN6+Z2JL9gg+Q7ybZ9VXy5bT+3nN3P7k0wpoP8+NuDUJQ/frbN7Sh+wz6dOsnzy/IJDwnixolpbkcxxm+k9Ijk5kn9mJ9VyM79h9yO4xesKHSC8pp6FuYUcs24ZOKjw92OY4xf+dfpAwgNFv7fJ1vdjuIXrCh0gldWFVDb0MxtU9LdjmKM30mMieD2KRksWltEXkmV23F8nhUFL6ttaOL5pflMHdSLoX2OuRfPGNMB/uVbGYQGB/Hc0ny3o/g8KwpetjBnD6VVddx57gC3oxjjt+Kjw7l6bF8W5hRSXlPvdhyfZkXBi5qalblLtjM6pTuTB8S7HccYv3b7lAxqG5p5ddVut6P4NCsKXvTBhr3kl9Vw57kDbJ4jY7xsWFIskwfE8/fl+TaD6hmwouAlqsozX2wno1cUF4/o43YcYwLC7VMyKK6o5cONtt7C6bKi4CXLtpexfk8Fc6b1JzjIWgnGdIbzhybSLz6S5+2C82mzouAlzy/dSXxUmK3BbEwnCg4Sbj0nnexdB1m7u9ztOD7JioIX7Co7xKebS7h5UhoRocFuxzEmoMzMTCE6PITnl+50O4pPsqLgBS8s20WwCN85u5/bUYwJODERoVw3Ppn31u/lwCEbnnqqrCh0sOq6Rl7P2s1lo5NIjI1wO44xAenms/tR39TM61k2PPVUWVHoYAuyC6mqa+T2KRluRzEmYA3uHcPE9J68sqqA5mZ1O45PsaLQgZqblReW5TM2NY6xqXFuxzEmoN18dhq7ymr4Km+/21F8ihWFDrRkm2dltdtt4jtjXHfJyD7ER4Xx0opdbkfxKVYUOtDfluWTGGMrqxnTFYSHBDMzM5VPN5dQXHHY7Tg+w4pCB9lRWs3nW0r5jq2sZkyXcdPENJpVmWfzIbWbfXp1kL8v30VYsK2sZkxXkhYfybRBCcz7uoD6RpsPqT2sKHSAqtoGXs/azeVjkkiIsZXVjOlKbpuSzr7KOt7ILnQ7ik+wotAB3sgu5FB9E7dPtmGoxnQ15w1OYGxqHE8tzqOuscntOF2eFYUzdGQY6oR+PRiV0t3tOMaYo4gI9184mD3lh5mfZa2Fk/FaURCRVBFZLCK5IrJRRO5xto8VkRUiskZEskRkYqvnPCwieSKyRUQu9la2jvTF1lLyy2q4bXK621GMMccxdVAvMvv14KnP8qhtsNbCiXizpdAIPKCqw4CzgbtFZDjwO+CXqjoW+IXzPc6+2cAI4BLgaRHp8rPJPb8sn96x4Vwy0tZMMKarOtJa2FtZy7xVBW7H6dK8VhRUtVhVc5zHVUAukAwocGQF++5AkfP4KmCeqtap6k4gD5hIF5ZXUs2SraV8Z1I/QoOtJ86YruycAfFMyujJU59vt9bCCXTKJ5mIpAPjgJXAvcDvRWQ38N/Aw85hyUDrwcSFzrajzzXH6XbKKi0t9Wbsk/r78nzPMNRJNgzVmK7uSGuhtKqO+TZR3nF5vSiISDSwALhXVSuBu4D7VDUVuA949sihbTz9mJmsVHWuqmaqamZCQoK3Yp9UZW0Db2QXcsWYvvSKtmGoxviCSf3jGZPSnb8v34WqTZTXFq8WBREJxVMQXlbVhc7mW4Ejj1/nn11EhUBqq6en8M+upS7npRW7qKlvsnmOjPExt5yTTl5JNct3lLkdpUvy5ugjwdMKyFXVJ1rtKgLOdR6fD2xzHi8CZotIuIhkAIOAVd7KdyYO1TXy1y93ct6QBEYm2zBUY3zJ5aOTiIsM5cXlNlFeW0K8eO4pwC3AehFZ42z7GfB94EkRCQFqgTkAqrpRROYDm/CMXLpbVbvk1aCXV+7iwKF6fnT+ILejGGNOUURoMDdkpvLXr3ZSXHGYpO7d3I7UpXitKKjqV7R9nQBgwnGe8xjwmLcydYTD9U3MXbKDqYN6MaFfD7fjGGNOw3fO7sfcL3fw6soC7r9oiNtxuhQbR3mKXllVwP7qen78bWslGOOrUntGMn1IIq+s2m0T5R3FisIpOHConj99to1z+sdzVnpPt+MYY87ALef0Y391HR9u3Ot2lC7FisIpeOzdXKpqG3n0yhFuRzHGnKFzByXQLz7SVmY7ihWFdlqWt58FOYXMmdafIX1i3I5jjDlDQUHCrMxUVu48QP7+Q27H6TKsKLRDbUMTj7y1gX7xkXYtwRg/cv2EFIIEu8O5FSsK7fDCsnx27j/Er68aSURol5+jzxjTTr1jI5g+JJE3sgtpbLILzmBF4aTKa+p5anEe5w1JYNpg96bVMMZ4xw1npVJSVcfnW9ydS62rsKJwEk8tzqOqrpGHLh3qdhRjjBdMH5pIr+hwXrMuJMCKwgntPlDDC8t2cd34FIb2iT35E4wxPic0OIjrJiTz2eYSSqpq3Y7jOisKJ/DEx1sRgfsvHOx2FGOMF83KTKWpWVmYs8ftKK6zonAcG/ZU8NaaPdw+JYO+cTY3ijH+bEBCNGel92B+1u6An1LbisJx/PaDzXTvFspd5w1wO4oxphNcOz6FHaWH2FRc6XYUV1lRaMOSraV8uW0/P5w+kO7dQt2OY4zpBBeP6ENwkPDe+mK3o7jKisJRmpuVx9/fTEqPbtxyTj+34xhjOknPqDAmD4jn3XXFAd2FZEXhKG+v3cOm4kp+cvEQwkPsRjVjAsmMUUnkl9UEdBeSFYVWahua+O8PtzIyOZYrRvd1O44xppMd6UJ6d13gdiFZUWjlxeW72FN+mIcvHUZQ0PHWBzLG+KuWLqT1gduFZEXBUVHTwJ8W53Hu4ASmDOzldhxjjEsuG5XErrIaNhYFZheSFQXH05/nUVnbYNNZGBPgWrqQAnQUkhUFYE/5YZ5fls8145IZlmTTWRgTyHo4XUjvBWgXkhUF4A8fbQHgAVvA2xgDXDG6L7vKalizu9ztKJ0u4IvCpqJK3ly9h9snp5Ns01kYY4BLRvUhPCQoIOdC8lpREJFUEVksIrkislFE7nG2vyYia5yvfBFZ42xPF5HDrfY9461srT3+wWZiI0L51/MGdsaPM8b4gNiIUC4e0YdFa4uoa2xyO06nCvHiuRuBB1Q1R0RigGwR+VhVbzhygIj8Aaho9ZztqjrWi5m+4att+1mytZRHZgyje6RNZ2GM+afrJqSwaG0Rn+WWcOmoJLfjdBqvtRRUtVhVc5zHVUAukHxkv4gIMAt41VsZTpKPxz/IJTnOprMwxhzrWwN7kRgTzoKcQrejdKpOuaYgIunAOGBlq81TgX2quq3VtgwRWS0iX4jI1OOca46IZIlIVmnp6S+f92luCRv2VHLvBYNs3WVjzDGCg4RrxiXz+ZZS9lfXuR2n03i9KIhINLAAuFdVW98NciPfbCUUA2mqOg64H3hFRI4ZH6qqc1U1U1UzExJOb81kVeV/P9tGWs9Irh6XfPInGGMC0nUTUmhsVt5eU+R2lE7j1aIgIqF4CsLLqrqw1fYQ4FrgtSPbVLVOVcucx9nAdsArS559vrWUdYUV/Ot5AwgNDvgBWMaY4xjcO4ZRyd1ZkB04XUjeHH0kwLNArqo+cdTuC4DNqlrY6vgEEQl2HvcHBgE7vJHtqc/ySI7rxrXjU7xxemOMH7lufDKbiivZsrfK7Sidwpt/Jk8BbgHObzXMdIazbzbHXmCeBqwTkbXAG8Cdqnqgo0Ot3V1O1q6D3PGtDMJCrJVgjDmxy0b3JUjgnbWB0YXktSGpqvoV0OZUo6p6WxvbFuDpavKq55fuJDo8hFmZ1kowxpxcQkw4kwf04p11RTxw0WA8nSD+K6D+VN5XWcu764uZmZlCTITdl2CMaZ8rxnhmTt2wx/9nTg2oojD/6900Niu3TU53O4oxxodcPKIPocHCO+v8vwspoIrCuj0VDEyIpl98lNtRjDE+JC4yjGmDEvjH2iKam/175tSAKgrbS6sZkBDtdgxjjA+6YkxfiipqySk46HYUrwqYotDQ1ExBWQ39E6yVYIw5dRcM7014SJDfj0IKmKJQcKCGxma1loIx5rREh4dw/tBE3l2/lyY/7kIKmKKwo/QQgLUUjDGn7YoxfdlfXcfKHWVuR/GagCkK20urAehvLQVjzGmaPiSRqLBgvx6FFDBFYUdpNb2iw+neze5PMMacnm5hwVw4vDfvb9hLfWOz23G8ImCKwvbSQwywriNjzBm6YkxfymsaWJq33+0oXhEwRWFHabV1HRljztjUQQnERoT47SikgCgKBw7Vc7CmwVoKxpgzFhYSxKUjk/ho0z5qG/xv/eaAKAo793suMttwVGNMR7h8TBLVdY18vuX0V3/sqgKiKJRWeZbSS4wNdzmJMcYfnNM/nvioML8chRQQReFgTQMAPSLDXE5ijPEHIcFBzBiVxKe5+6iua3Q7TocKkKJQD0BcpA1HNcZ0jGvHJ1Pb0Ox3S3UGRFGoqGkgLCSIbqHBbkcxxviJcWk9GJcWx/NLd/rVzKkBURQO1tTTIzLU71dMMsZ0rjumZJBfVsPiLSVuR+kwAVIUGux6gjGmw10ysg9J3SN49qudbkfpMAFRFMpr6u16gjGmw4UGB3Hr5HSWbS8jt9g/luoMiKJgLQVjjLfMPiuVbqHBPL/UP1oLAVEUPC0FKwrGmI4XFxnGdROSeWtNEfur69yOc8a8VhREJFVEFotIrohsFJF7nO2vicga5ytfRNa0es7DIpInIltE5OKOyKGqlNc00MO6j4wxXnL7lAzqG5t5eUWB21HOWIgXz90IPKCqOSISA2SLyMeqesORA0TkD0CF83g4MBsYAfQFPhGRwap6RpOLVNU10tis1n1kjPGaAQnRTB+SwIsrdnHnef0JD/Hd4e9eaymoarGq5jiPq4BcIPnIfvGMD50FvOpsugqYp6p1qroTyAMmnmmOCuduZrvQbIzxpju+lcH+6jreWVvsdpQz4s2WQgsRSQfGAStbbZ4K7FPVbc73ycCKVvsLaVVEWp1rDjDH+bZaRLa0J8Os355i6H/qBfjnxOnHF4iv2Rvs93jmfO53OPP0P2ta8/br7ne8HV4vCiISDSwA7lXV1mO2buSfrQSAtu4sO+Y2QVWdC8zt0JAnICJZqprZWT+vKwjE1+wN9ns8c4H6O3TzdXu1KIhIKJ6C8LKqLmy1PQS4FpjQ6vBCILXV9ymA/01BaIwxXZg3Rx8J8CyQq6pPHLX7AmCzqraeSWoRMFtEwkUkAxgErPJWPmOMMcfyZkthCnALsL7VsNOfqep7eEYZte46QlU3ish8YBOekUt3n+nIow7SaV1VXUggvmZvsN/jmQvU36Frr1tU/Wd2P2OMMWcmIO5oNsYY0z5WFIwxxrSwonAUEZnpTMvRLCLHHRImIpc403HkichDnZmxo4lITxH5WES2Of/2OM5x+SKy3pmiJKuzc3ZFJ3sfiMf/OvvXich4N3J2de34PZ4nIhWtpsj5hRs5O5KIPCciJSKy4Tj7XXnvWFE41gY8w2WXHO8AEQkGngIuBYYDNzrTdPiqh4BPVXUQ8Knz/fFMV9WxgTh2/GjtfB9cimck3SA8N13+X6eG9AGn8P/Tl857b6yq/qpTQ3rH34BLTrDflfeOFYWjqGquqp7sLumJQJ6q7lDVemAenmk6fNVVwAvO4xeAq92L4lPa8z64Cvi7eqwA4kQkqbODdnH+9v9Tu6jqEuDACQ5x5b1jReH0JAO7W33f5pQcPqS3qhaDZ84qIPE4xynwkYhkO9ONBLr2vA/87b3iDe39HZ0jImtF5H0RGdE50VzlynunU+Y+6mpE5BOgTxu7HlHVt9tzija2demxvSd6zadwmimqWiQiicDHIrLZ+WsnULXnfeBz7xUXtOd3lAP0U9VqEZkBvIWnW8WfufLeCciioKoXnOEpfG5KjhO9ZhHZJyJJqlrsNE/bXIVcVYucf0tE5E08zf5ALgrteR/43HvFBSf9HbWeN01V3xORp0Wkl6r61GR5p8iV9451H52er4FBIpIhImF47tBe5HKmM7EIuNV5fCtwTGtJRKKcdTEQkSjgIjwX5QNZe94Hi4DvOiNJzgYqjnTVmRYn/T2KSB9n6hxEZCKez66yTk/auVx57wRkS+FEROQa4I9AAvCuiKxR1YtFpC/wV1WdoaqNIvJD4EMgGHhOVTe6GPtMPQ7MF5HvAQXATIDWrxnoDbzp/H8ZAryiqh+4lLdLON77QETudPY/A7wHzMCzPkgNcLtbebuqdv4erwfuEpFG4DAwW318OgYReRU4D+glIoXAfwCh4O57x6a5MMYY08K6j4wxxrSwomCMMaaFFQVjjDEtrCgYY4xpYUXBGGNMCysKxrQiItUn2R8nIv/a6vu+IvKG83isc7ftqf7MR0Xk3049rTEdz4qCMacmDmgpCqpapKrXO9+OxTOu3BifZUXBmDaISLSIfCoiOc4aEkdm7XwcGODM6f97EUkXkQ3Onbi/Am5w9t1wdAvAOS7defyIs37AJ8CQVscMEJEPnEkHvxSRoZ33qo2xO5qNOZ5a4BpVrRSRXsAKEVmEZ62Jkao6FuDIh7yq1jsLv2Sq6g+dfY+2dWIRmYBnKodxeP4fzAGynd1zgTtVdZuITAKeBs73yis0pg1WFIxpmwC/EZFpQDOeKYt7d9C5pwJvqmoNgFNsEJFoYDLwujOdCEB4B/1MY9rFioIxbbsZz/xXE1S1QUTygYhTPEcj3+yibf38tuaXCQLKj7RCjHGDXVMwpm3dgRKnIEwH+jnbq4CY4zzn6H35wHgAZ33dDGf7EuAaEenmzDx7BbRMD71TRI5MSCgiMqbjXpIxJ2dFwZi2vQxkikgWnlbDZgBVLQOWOheNf3/UcxYDw49caAYWAD1FZA1wF7DVOUcO8Bqwxjnmy1bnuBn4noisBTYSAMtSmq7FZkk1xhjTwloKxhhjWlhRMMYY08KKgjHGmBZWFIwxxrSwomCMMaaFFQVjjDEtrCgYY4xp8f8Bbz9OB0ckyJkAAAAASUVORK5CYII=\n",
"text/plain": [
"