{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Selecting and Reindexing of Area of Interest" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### _The elegant way of dealing with large output_" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "_Before working through this script, it is helpful to have had a look into [Triangular Meshes and Basic Plotting](tripcolor.ipynb) to get a basic understanding of plotting on a triangular basis._\n", "\n", "Three-dimensional global ICON output for example requires many times more memory than two-dimensional output. The handling of such large amounts of data can very quickly lead to the exhaustion of the given memory. If you are sitting right next to a supercomputer, you are tempted to just request more RAM and go for it. However, there are elegant solutions besides the powerful one and since a lot of RAM also means a lot of electricity and a lot of coolant, the motivation for an elegant way is also to save limited resources.\n", "\n", "In many cases, we rarely look at the complete global output, but rather at a specific, selected area. It is therefore advisable to cut out only this area. For this purpose it is very helpful to reduce the grid information from the global grid file to the area of interest but in such a way that the indexing makes sense starting at 0 and counting up continuously. The advantage is that we generate a new local grid-file, which looks like the global grid-file but is much smaller in terms of storage capacity and therefore easier and faster to handle.\n", "\n", "So let's do something for the climate 🌍 and first of all load the necessary libraries and the global grid-file:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import xarray as xr\n", "import numpy as np" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Importing the Grid-File" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
       "Dimensions:                         (cell: 20971520, nv: 3, vertex: 10485762,\n",
       "                                     ne: 6, edge: 31457280, no: 4, nc: 2,\n",
       "                                     max_stored_decompositions: 4, two_grf: 2,\n",
       "                                     cell_grf: 14, max_chdom: 1, edge_grf: 24,\n",
       "                                     vert_grf: 13)\n",
       "Coordinates:\n",
       "    clon                            (cell) float64 ...\n",
       "    clat                            (cell) float64 ...\n",
       "    vlon                            (vertex) float64 ...\n",
       "    vlat                            (vertex) float64 ...\n",
       "    elon                            (edge) float64 ...\n",
       "    elat                            (edge) float64 ...\n",
       "Dimensions without coordinates: cell, nv, vertex, ne, edge, no, nc,\n",
       "                                max_stored_decompositions, two_grf, cell_grf,\n",
       "                                max_chdom, edge_grf, vert_grf\n",
       "Data variables: (12/91)\n",
       "    clon_vertices                   (cell, nv) float64 ...\n",
       "    clat_vertices                   (cell, nv) float64 ...\n",
       "    vlon_vertices                   (vertex, ne) float64 ...\n",
       "    vlat_vertices                   (vertex, ne) float64 ...\n",
       "    elon_vertices                   (edge, no) float64 ...\n",
       "    elat_vertices                   (edge, no) float64 ...\n",
       "    ...                              ...\n",
       "    edge_dual_normal_cartesian_x    (edge) float64 ...\n",
       "    edge_dual_normal_cartesian_y    (edge) float64 ...\n",
       "    edge_dual_normal_cartesian_z    (edge) float64 ...\n",
       "    cell_circumcenter_cartesian_x   (cell) float64 ...\n",
       "    cell_circumcenter_cartesian_y   (cell) float64 ...\n",
       "    cell_circumcenter_cartesian_z   (cell) float64 ...\n",
       "Attributes: (12/43)\n",
       "    title:                    ICON grid description\n",
       "    institution:              Max Planck Institute for Meteorology/Deutscher ...\n",
       "    source:                   git@git.mpimet.mpg.de:GridGenerator.git\n",
       "    revision:                 d00fcac1f61fa16c686bfe51d1d8eddd09296cb5\n",
       "    date:                     20180529 at 222250\n",
       "    user_name:                Rene Redler (m300083)\n",
       "    ...                       ...\n",
       "    topography:               modified SRTM30\n",
       "    subcentre:                1\n",
       "    number_of_grid_used:      15\n",
       "    history:                  Thu Aug 16 11:05:44 2018: ncatted -O -a ICON_gr...\n",
       "    ICON_grid_file_uri:       http://icon-downloads.mpimet.mpg.de/grids/publi...\n",
       "    NCO:                      netCDF Operators version 4.7.5 (Homepage = http...
" ], "text/plain": [ "\n", "Dimensions: (cell: 20971520, nv: 3, vertex: 10485762,\n", " ne: 6, edge: 31457280, no: 4, nc: 2,\n", " max_stored_decompositions: 4, two_grf: 2,\n", " cell_grf: 14, max_chdom: 1, edge_grf: 24,\n", " vert_grf: 13)\n", "Coordinates:\n", " clon (cell) float64 ...\n", " clat (cell) float64 ...\n", " vlon (vertex) float64 ...\n", " vlat (vertex) float64 ...\n", " elon (edge) float64 ...\n", " elat (edge) float64 ...\n", "Dimensions without coordinates: cell, nv, vertex, ne, edge, no, nc,\n", " max_stored_decompositions, two_grf, cell_grf,\n", " max_chdom, edge_grf, vert_grf\n", "Data variables: (12/91)\n", " clon_vertices (cell, nv) float64 ...\n", " clat_vertices (cell, nv) float64 ...\n", " vlon_vertices (vertex, ne) float64 ...\n", " vlat_vertices (vertex, ne) float64 ...\n", " elon_vertices (edge, no) float64 ...\n", " elat_vertices (edge, no) float64 ...\n", " ... ...\n", " edge_dual_normal_cartesian_x (edge) float64 ...\n", " edge_dual_normal_cartesian_y (edge) float64 ...\n", " edge_dual_normal_cartesian_z (edge) float64 ...\n", " cell_circumcenter_cartesian_x (cell) float64 ...\n", " cell_circumcenter_cartesian_y (cell) float64 ...\n", " cell_circumcenter_cartesian_z (cell) float64 ...\n", "Attributes: (12/43)\n", " title: ICON grid description\n", " institution: Max Planck Institute for Meteorology/Deutscher ...\n", " source: git@git.mpimet.mpg.de:GridGenerator.git\n", " revision: d00fcac1f61fa16c686bfe51d1d8eddd09296cb5\n", " date: 20180529 at 222250\n", " user_name: Rene Redler (m300083)\n", " ... ...\n", " topography: modified SRTM30\n", " subcentre: 1\n", " number_of_grid_used: 15\n", " history: Thu Aug 16 11:05:44 2018: ncatted -O -a ICON_gr...\n", " ICON_grid_file_uri: http://icon-downloads.mpimet.mpg.de/grids/publi...\n", " NCO: netCDF Operators version 4.7.5 (Homepage = http..." ] }, "execution_count": 2, "metadata": {}, "output_type": "execute_result" } ], "source": [ "grid = xr.open_dataset(\n", " \"/work/mh0287/k203123/Dyamond++/icon-aes-dyw2/experiments/dpp0029/icon_grid_0015_R02B09_G.nc\"\n", ")\n", "grid" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since Max Planck was an important person in science, we will investigate the area around his birthplace. According to [this wiki](https://kiel-wiki.de/Max_Planck#Ehrungen_durch_die_Universit.C3.A4t) (sorry, only available in german) of the city of Kiel, he was born in a house (which unfortunately no longer exists) under these coordinates: 54.32 N, 10.13 E. We take this point as center and choose a 0.5°x0.5° rectangular window around it." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "max_plancks_birthplace_x, max_plancks_birthplace_y = np.array([10.13, 54.32])\n", "\n", "left_bound = 10.13 - 0.25\n", "right_bound = 10.13 + 0.25\n", "top_bound = 54.32 + 0.25\n", "bottom_bound = 54.32 - 0.25" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's see if we're right:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/sw/spack-levante/mambaforge-4.11.0-0-Linux-x86_64-sobz6z/lib/python3.9/site-packages/numpy/core/_asarray.py:102: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.\n", " return array(a, dtype, copy=False, order=order)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAATcAAAEVCAYAAACFRqIPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABd0klEQVR4nO2deXhNRx/HP5OIxBZEEmuESCJB7aJ2oZZSW1VRSku9ba1V6m0tRasoWrWV1lJVaqmtliIIoqutllgiCUHEFlsSWWT5vX/cK2+QSEKSm1zzeZ7z5Nw5s3zPyT3fOzNnZo4SETQajcbcsDC1AI1Go8kOtLlpNBqzRJubRqMxS7S5aTQas0Sbm0ajMUu0uWk0GrNEm9tzglJKlFKueakMpdQEpdTyJxw/qZRqnsk89yql3skOPZrchTa3bEIpFaKUuq+Usn8k/KjRBCpkcXkVjPlGGbcQpdTHWVnGs6CUyvIBlSJSVUT2PqFMbUbPMdrcspfzQM8HH5RSLwAFsrnMYiJS2Fjup0qpttlcXq5EKZXP1Bo0pkWbW/byE9Anxee+wLKUEZRS7ZVS/yqlIpRSl5RSE1Ic666UOqeUsjV+flkpdVUp5ZBewSLyF3ASqPbosXTKfFAD7KuUuqiUCldKjUlx3FIpNVopFayUilRKHVZKOaVSRmNj3t6pHHvLeF6RSqnzSqleTzgVG6XUamPcI0qpGinyCVFKvWTcn6CUWquUWq6UigDeA0YD3Y012WMp8nRWSv1hzNPnQe06xbn/RykVppS6opQakZYwpdQvxv/HXaWUn1KqaopjBZRSXymlLhiP/66UKmA89qJS6k+l1B2l1LHMNq01GURE9JYNGxACvAQEAJ6AJXAJcAYEqGCM1xx4AcMPTXXgGtA5RT4rgKVACSAMeCWN8ioY880HKKAREA20NB4XwDW9MlPksxBDLbMGEAd4Go9/BJwAKhvLqQGUSFkG0MZ4rl6p6CwERACVjZ9LA1XTOKcJQDzwGmAFjMRQG7ZKeY0fidvZeF4FjGHLH8lzLxAMuBvj7AWmPnLuK406XwBuPFLG8hR59QOKANbAN8DRFMfmGfMua/zfNzTGKwvcBNoZdbYyfnYw9XfW3DaTCzDXjf+b21hgCtAW2Gk0n2RzSyXdN8DMFJ+LAReNhvLdE8p7cGPeAW4Dp4GhKY4nm9uTykyRT7kUxw8APYz7AUCnNPIR4BPgAvBCGnEKGTV2BQqkcw0nAH+n+GwBXAGapLzGKeL6pZI+NXMbm+LzQGD7I+fukeL4NGBxWvk98n8SoKhRZwxQI5V4/wV+eiRsB9DX1N9Zc9t0szT7+Ql4A3iLR5qkAEqp+kqpPUqpG0qpuxiaU8kPIUTkDvALhublVxkoz15EiouIp4jMTi1CemUauZpiPxoobNx3wlDzSYsPgDUiciK1gyJyD+huLPOKUmqrUsrjCfldSpE2CQgFyqQXNx3SOrfU8rmQWnnG5vlUY/M8AoPRguE62gM2pH6dnIFuxibpHaXUHaAxhhqsJgvR5pbNiMgFDE2pdsD6VKL8DGwCnESkKLAAQ3MPAKVUTQzNn5VAqmb1FDyxzHS4BFR6wvFuQGel1AdpRRCRHSLSCsMNfQZDEzgtkvvzlFIWQDkMzfNUs07nc0ZJ2YdYPo3y3gA6YaidF8VQ6wPDdQwHYkn9Ol3CUHMrlmIrJCJTn1KrJg20ueUM/YEWxlrLoxQBbolIrFLKC8NNA4BSygZYjqFj/G2grFJqYBboSbPMDLAI+Fwp5aYMVFdKlUhxPAxoCQxNTatSqqRSqqNSqhCGvrwoIPEJ5dVRSr1qfPr5gTHN3xnUeg2oYDTFzDBOKVXQ+IDgbWB1KnGKGLXcBAoCkx8cMNYwlwBfK6XKGGt5DZRS1hj+nx2UUm2M4TZKqeZKqXKZ1KhJB21uOYCIBIvIoTQODwQ+U0pFAp8Ca1IcmwKEish8EYkDegOTlFJuzyjpSWWmx9fG+D4YHgws5pHhLSJyEYPB/Vc9PmDWAhiBwQRvAc2MetLiVwzN2NvAm8CrIhKfQa2/GP/eVEodyWAagH1AELAbmCEiPqnEWYahyXoZOMXjhjsSQz/pQQzn+SVgISKXMNT4RmN4WHEJw0MafS9mMcrYoanRPPcow8DqB09jE0wsR/OM6F8LjUZjlmhz02g0Zolulmo0GrNE19w0Go1Zos1No9GYJWZjbkqpYUopf2VY4+sDY1gZpZSvUupXpVRhY9gEpdRlZVh66MFWLKe1PXK8uXFy9QM9n6Y41sM4YfyDFGEhSqkTKeJn1eDelJqWKKWuK6X8U4TZKaV2KqUCjX+Lpzg2XSl1SCnVzPi5glIq5pHr3Ce1srJT1yNpU163QynCc/R7ksY5dDN+P5KUUnUfiZ8j19bsMPX8r6zYMExN8scwmDIfsAtwA6YCVYEOwHvGuBOAkabW9kic5sCWNNJvxDDxehVQ2BgWgmGaVXbqbgrUBvxThE0DPjbufwx8adz3AKYbz3GNMaxCyrSm0JVK2lSvW05/T9I4B08MixHsBeqmCM+xa2tum7nU3DwxTLCOFsP4pH1AFwymkGTcMjq9KKe0ZZQHuoUcPAcR8cMw+DQlnYAfjfs/YliBA/5/nbNdYyZ1ZZQc/Z6kdg4iclpEAp6gLUf//+aAuZibP9BUKVVCKVUQwzxOJ2Au8B2GSdopV2QdnqI6v8dE2h6lgTKs7bVNpVgXDMN81EPAIRGJTBG+J8U5DM8++Q9RUkSuABj/Ohr3T2KoWfwOzE8Rv9IjTacmOakrFQTwUYY16P6TIjw3fE9SJRdc2zyLWaxWKiKnlVJfYlhSKAo4BiSIYdJ601SSzBSRGabU9ki0I4CziEQppdphaIq6GdP/yP9rJSnxFpHwbBOeSURkSCrBwSJSM6e1PIFGIhKmlHIEdiqlzoiIX274njyJPHJtcx3mUnNDRBaLSG0RaYqhyh9oak0PSE+biESISJRx/zfASj3y7oVcwjWlVGkA49/rJtbzgAzpEpEw49/rwAbAK8cUanIcszE3468xSqnywKsYlgjKFaSnTSlVSimljPteGP4vN3NaZwbYhGGpdIx/fzWhlpSkq0spVUgpVeTBPtAaQ5eBxlwx9RONrNqA/RhWZziGcWntNOJNwLCSw9EUW4Wc1oahf+fBk7nBGN53cAzD6hIN08kvBMOKEw/0L8sGzSsxrHobj2GByP4YljrfjaHmuRuwe0L6ChhWo015nYfmpC4Mi0z+Ztx3MV7fY8ZrPSadcrLte5LGOXQx7sdhWKppR05fW3Pb9PQrjUZjlphNs1Sj0WhSos1No9GYJdrcNBqNWaLNTaPRmCXa3DQajVmSp2cotG3bVsLDMzZI//79++TPnz+bFWUNWaU1IiKCGzdukJCQQGxsLNbW1jg7O1OgQIH0E2eClHpjYmKIjo4mOjo6ed/CwoICBQpQsGDB5M3a2jpLNTwgOjqa4OBgHB0dKVmy5BO1PoqIcOSI4T0yjo6OODk9PEsuNDSUu3fv4u7ujpWVVdaLT4XMfBdiY2MJDAzkhRdeyFINhw8f3iEibbM00xwgT5tbeHg4hw6l9VKph9m7dy/NmzfPXkFZRFZoTUxMpFixYsyZMwcXFxecnJyoWLEiMTExhIaGcvHiRS5duoSI4OLiQtOmTTGOI840O3bswM/Pj6VLl1K4cGFefPFFatasSa1atahZsyaOjmlN9cweLl68SOvWralQoQKtWrXC3t6eatWq4ebmhp+f3xOv7aRJk/Dx8WHlypWULVv2oWMiwhdffMGPP/7Ivn37KFMmrXdDZx2Z+S4cOXKEVq1aZfieyCi5dLZM+mRw0GEI/x80eiiNQY7tMpo2xQBLXwyjyQunyDMacEwRLyotXXXq1JGMsmfPngzHNTXPqvXmzZvy22+/iZ2dnYiIJCUlydy5c6VkyZJibW0tLi4u0qxZM+ndu7f06dNHSpYsKZs2bcpUGdHR0bJ8+XLx9vaWb775RoYPHy7+/v7PpDsriYyMlNmzZ8s777wjnTp1EgyT5uW7776TM2fOPDFdqVKl5NSpU2nGeffdd+Xzzz/PDtmP8aTvwtWrV2XXrl3y008/ydChQ6VcuXKybNmyLNeQ8r6VPOALD7bM1NxSm6g9UzI2sTi1tEOBIRhGjvfG8NZzMLytewTw30xoM3uuXr1KUlISpUuXfqiGdfnyZZYvX87x48cJCgoiMDCQxMRE3NzcGD9+PLdu3aJ///5cvHiR3bt34+npiYXFw12tEydOZMWKFbzyyivp1t6OHj3K4sWL+fnnn6lbty7vv/8+9vb2DBs2LFvO+2kpXLgwQ4YY5pt/8803/PHHH1haWhIZGYmHhwehoaHJNTNfX18+/vhjXnnlFT788ENatmzJ/v378fT0TDXvhg0bsm7dOsMo+Kes7T4rq1ev5t1336VWrVo4OjpSt25dVq9eTcOGDXNaSq71BVM+UEhrDa0lQHellJ1JVOVCwsPDKV26NGXLlsXJyYkePXrw4Ycf0q5dO1544QXOnz9PmzZtmDlzJgEBAdy5c4ddu3ZRokQJatWqRcWKFfnzzz+pWrXqY8YG8NFHH3Hp0iWcnZ0ZPHgwu3btIj4+nqSkJM6ePcvPP//MiBEjqF27Nh07dqREiRIcOXKEHTt20K1bN5Pd4Bll3rx5TJ8+natXr1KnTh2SkpIeanIePHiQoKAgduzYgaOjI1u3bqVWrVpp5te+fXvOnz9Pp06dCAw0zfoMs2fPZuXKlSxdupS6detiaWnJ5cuXuX37tkn0ZCFZ5wvpVe2MVcDzGJblOQz8J0VVMQQ4biy4eEbTGsOdAT9gM1AkRZ4jMbwFfWJ61U9zbpYmJSXJpUuXZMuWLTJw4EBRSklwcLAEBgbKTz/9JNOmTZO1a9dKRETEQ2n//fdf8fLyksKFC0u7du1k+/btGSozKSlJTp48KZMnT5b69etLsWLFxNbWVpydnaVr164yZcoU8fX1lYSEhFT15mZWrFghZcuWlVGjRqWqNSYmRsaMGSMODg4CyFdffZVunjExMTJ16lQpUaKE7Ny5MxtUG0jr2jZv3lzc3d3Fzs5O3n33XRk2bJi0b99ebG1tZcSIERIdHZ1lGki7WZorfeHBlqG5pUqpMpJiHSwM1cYADFVFAT4HSotIv4ykFcNKpKmVMwHDmmeLMLTFqwNhIlI4tfjVq1eX2bMz9vqAqKgoChdONZtcRVJSEnfu3OHSpUsAyU8XHR0dSUxMJDo6msjISCwsLChWrBiFCxd+qObk7++PUgoHBwesra2xtbV9qppVQoJhybl8+dLvucjN1/bBE9uIiAiSkpIoWbLkE7XevXuXc+fOPbHmlpJbt24RHh6Ou7t7Vkl+iLSubWJiIrGxsRQsWPCh/29CQgIXLlwgMTERV1fXVGvqmcXb2/uwiNR9NDy3+sIDMtTnJinWwVJKbQC8UgpRSi0EtmQ0LQZnflJ5d5RSPwMDnxQvf/78GX6SlFeelq5cuZKwsDBWr16Np6cnN27c4OzZs4SGhuLk5ETVqlVp3bo1t2/fZs6cOQQHB/Pyyy/TqVMn2rVrR0xMDP7+/uzbtw9/f3+Cg4MZPXo0/fv3z7bhF7nx2sbHx/Phhx+yYcMGmjdvntzkDggIeEjr2bNn2bp1K//88w/Hjx/n/PnzNG7cmGHDhj3RGESEo0ePMm/ePBITE/nPf/6TZtxn4WmubWJiIr179+bYsWPMmzcvW3RB7vWFlAnSa5IW4v/Vw0LAn0BbDI78IM5wYFVG0z6hrAkYX8oB2GOousamFd9cm6U+Pj6yePFiWbJkiWzatEnOnDkjcXFxqcYNDQ2V+fPnS6NGjaRv377J4cHBwRIQECAHDhyQdu3aiZOTkyxYsCDNfJ6F3HZtw8PDxdvbW15++WW5c+dOcnhSUpL4+vrK7t27Zfjw4eLm5ialS5eWAQMGyLJly+TYsWNPvD5xcXHy+++/y6RJk+SFF14QZ2dnmThxokRGRmbbuTzttb1586Y4ODjIyZMnn1kDqTRLc7MvJKdLN0Ia62ABP2F4lHscw2KBpY3hz7qG1sgUn782+O/zZW5Po/XAgQPi4eGR/Ll9+/YCSNmyZWXQoEGyZs0aadOmjTg7O8vChQvl/v37JtWbXZw/f15cXFxk1KhRD/UPLlu2TKytreXrr78WLy8vmThxohw+fFgSExNl1apV0q5dO6lUqVLykBFLS0upU6eOvPnmmzJo0CBp0aKFFC5cWGrVqiUffPCB7NmzRxITE7P9fJ7l2n799dfSrl27Z9aQhrnlWl9IjpdehNy8aXMzkJSUJL1795bhw4cnh50+fVpq1KghVapUkXfeeUfs7OxkwIABMnfuXHnppZekYsWKsmTJEomPj89xvdnJ6dOnxd7eXjZu3CiRkZESHh4umzZtEkdHRzl27Jjs2rXrofirVq0SNzc3Wb16tXTu3FkKFiwonp6eySaXcuvevbv88ssvOWJqD3iWaxsXFycVK1aUv/7665k0pGZueWHL0zMUnlfu3btHQEAAly5dIiQkhA0bNhAbG8v8+f9/OZKHhwf//PMPY8eO5eeff6Zbt274+vqycOFCwPCgol+/fhw4cOChdHkdDw8PfvvtNzp06EBkZCTW1tYUKVKEtWvXUr16dfbu3ftQ/L179+Lg4EBAQAAbN26kffv2bNny/26iuLg4Lly4wLlz5wgKCmLatGl8+eWXrF27Fmdn5xw+u8yRP39+BgwYwNKlS3nxxRdNLSfH0eaWx/jwww9ZtGgRFStWpHz58jg5OTFw4EA6der02AMDa2trpk+fTrdu3di8eTN16tShfPny7N69m+joaPLly0fr1q1NdCbZR7169bh69WqG4n7zzTcsXLiQq1evsnjxYl599dWHjltbW+Pu7p78NHTQoEF8/fXXNGjQgLVr15pi0GymaNq0abY97MjtaHPLY1hZWdGxY0eWL1+efmQjXl5eeHn9/0VPBw4coG3btty7d48ePXpw6NChLJ9snVewtrZm8ODBGY6vlGLEiBF4enrSuXNnpk2bxltvvZV9Ap+R6OhowsLCWLhwIQMGDDC1nBxFL3mUxxgzZgxHjhzh+++/f+o8LC0tuX37Nvfv3+f+/fu88847zJ07l/j4+CxUat60a9eOvXv3MmnSJD766CMSExNNLSlVWrVqxddff81vv/1maik5jja3PIatrS2rVq1i7NixT51HnTp1EBESEhL4559/KFq0KJ999hlRUVFZqNT8qVKlCv/88w+HDx+mY8eOREREmFpSMhERERw5coRvvvmGcePG5eraZXahzS0Pcvv2beztn24VmsjISLp27UqZMmWwtramQ4cOWFtbs337dooXL57FSs2fEiVKsGPHDpydnWnQoAHBwcEm0yIirF69msaNG1OmTBn69u3LmTNnWLx4MZ06dTKZLlOh+9zyIHPmzKFr166IZH5VioIFC2Jra4ulpSUHDhygdu3a2aTy+cHKyopvv/2Wb7/9lkaNGrFq1apsn7ERGRnJhg0bKFy4MIUKFSIiIoIlS5Zw9epVJkyYQNu2bbNtRkpeQZtbHqRPnz4MHz6ctWvX8vHHH/Pmm29meA6hpaUlP/zwA8uXL6d169Z06dKFzp07U6hQIe7cuUN0dDRFihShefPmFClSJJvPxLwYOHAglStXpnv37ixYsIAuXbpkW1l//fUXffv2pXPnzkRFRVGoUCEaNWpEvXr1sLCw4MSJE3h4eOTaOb85gqkH2j3L9jwP4n0wlcjLy0vatWv3VFOAwsPD5fPPP5eWLVtK06ZNpWPHjtK9e3d56aWXxNHRMcMrS5jbtX1WDh06JA4ODnL16tVnzistvffu3ZNKlSrJ1q1bRUQkLCxMnJycpHHjxvLSSy9JjRo1pGjRorJu3bpn1kAeHcRrcgHPsj3P5vaA+/fvS7169eT777/PsvITEhIkX758GZ6iZa7X9lkYOXKk9OvX75nzeZLeXbt2SdmyZWXmzJnSvHlzGTVqlERFRcnmzZtlwYIF0qhRIxkzZswza8ir5qYfKORxPvvsMyIiIrK07+zmzZsULVo0x16CYo6MGzeObdu2ceDAgWwro2XLlg+tDDN58mTmzJlDhw4d2Lt3L926deOjjz7KtvJzO7rPLY+jlCJfvnxZ2j929+5dbGxsEDHdMtp5lRs3bnD8+HHKly/PlClTGDJkCH/99VeWrKuWGl26dHmob2/EiBEEBQWxZcsW6tev/1w/VNA1tzzOxIkTGTJkCN7e3pw+fTpL8nR1daVEiRJs2rQpS/J7Xrh06RJubm5MnDiRJk2asGjRIgIDA9m5c2eOabCysmLRokX4+Pjg6+tL5cqVWbJkSfLio88T2tzyOEop3n33XaZMmULLli3x9/fPkjw///xzJk+enAUKnx8WL17Mm2++iZ+fH5cuXaJz587cvn2bX375Jce1VK9enU2bNrFy5UqWLVuGp6cnkyZNIiwsLMe1mAptbmZCnz59mDFjBq1atWLlypXcuXPnmfKrV68eJ06cMDx10mSIQoUKce3aNUQEKyur5CbitGnTTKapYcOG7Nmzhx9//JHLly9Tr149YmNjTaYnJ9HmZka88cYbLFq0iBUrVuDu7s7WrVufKp+kpCS6dOnCoEGDdJ9bJnj33XcJCgpi8ODByXNNK1WqhJ2daV/kppSiYcOGzJkzhxs3bmBpaWlSPTmFNjcz48F6ZLNmzSKjL895lKtXr3Ly5EnGjBmTxerMG1tbW/bu3cupU6d44403ct1CBNeuXcPCwoKDBw+yc+dOPvnkE6pVq2ZqWdmGNjczxd3dnQsXLjxV2jJlytC/f386d+6sm6WZxNbWlm3btnHlypXkhUFzC7a2tsTFxdGoUSMGDx7M1KlTc+1qJlmBNjczJC4ujpEjR9KjR4+nzuOrr74iNjaW9u3b4+vra9Y3QVZjY2PD1KlTmTNnjqmlPISFhQX58uXDxcWFu3fv0r17d3bt2mVqWdmGHudmhkydOhVbW1s+/fTTp87D0tKSffv2sWTJEkaMGMGZM2dwcXHB1dUVT09POnXqxIsvvqj75NKgatWqXL582dQyHqJQoULUq1ePSZMm0aJFC1PLyXa0uZkhu3btYsKECc88cNTa2pr333+f999/n3v37hEcHExgYCDHjh3jtdde49tvv30ul9LJCDY2NsTExOS6gdBvvPEGgwcPZseOHTg5OQGGWS4HDhx4bEZKQkICJ0+eNIXMLEGbWx7h7t27/Prrr4SFheHo6Ei5cuUoW7YsZcuWpWjRoskvY163bh0hISEPLSueFRQqVIjq1atTvXp12rVrx08//UShQoWYMGECFStWxMfHxyzfx/C0XLx4EScnp1xlbACDBw8mLi4OT09PatSoQc2aNfn2228BQ2095Zg8pRSVK1emSpUqppL7TGhzywP4+/vz8ssvU7duXVxdXTl9+jSXL1/m8uXLhIaGEhcXh4WFBZ6ennh7e7Nv375sXa7on3/+oWTJkuzbt4/9+/dTpUoVBg0ahIeHB4MGDaJAgQIcOHCA7du3c/LkScqVK8fLL7/Mp59++tzMV92+fXuufePUiBEjGDBgAIcPH+bEiRNMmTKFP//8E1dX12xdpimn0eaWizl9+jQ//fQTCxcuZPbs2fTs2TPVeLGxsVhbW+dYLaFIkSLcunWLwoULU6FCBRwdHfH392f27Nl89dVXREdHU7NmTYYPH06NGjW4cOECkydPpnv37qxZs4Z8+cz7axcREcHMmTP54YcfTC0lTWxtbfH29sbb29vUUrIN8/6W5THu3r3LggULKFmyJB4eHkRERPDGG2/wxx9/JL9aLjVsbGxyUCXUqlWLkiVLcvfuXQ4ePEjLli1p0qQJH330UaqrUDg5ObFhwwbat2/Pp59+atbTukSE9957j1atWtGsWbMMpfnzzz8JDw+nWbNmFC1aNJsVPj/ooSC5iKlTp7Jz506KFCnCqlWruHTpEjNmzHiisZkCCwsLFi9ezPfff8+mTZuIiYnh9ddfJyYmJs001tbWrFy5khUrVjBjxgyzHT/3999/c+DAAb755psMxV+3bh1dunRh9uzZlCxZkoMHD2avwOcIbW65hMjISAoUKMDu3bu5dOkS+/bty9XTZPLly8fNmzdJSEjAzc2N+Pj4dIeeODg48Pvvv7Nw4UImTZqUQ0pzluPHj+Pt7U2BAgXSjbtmzRoGDRrEjh07+OmnnyhQoACenp45oPL5QJubiYmPj8fLywtbW1tWr16NlZUViYmJ3L5929TSnkhAQAAA8+bN49q1a9y6dQtHR8cnpklISMDPz4/bt2/Tpk2bnJCZ4wQHB2doLunKlSsZNmwYPj4+1KxZk2XLltG1a9fn+50HWYzuc8sFFCxYkPz583P27FksLCxQSrFlyxZKly7NgAEDsm2hw2ehVatWrF27ltDQUBITE6lduzZvvvlmqnHv3LnDokWLmDNnDk5OTmzcuDHLh6qYmoSEBD755BN++eUXfH19nxj3wIEDfPDBB/j6+lK1alVEhCVLlrB06dKcEfu8YOp1zp9lM7d3KMTFxUlMTIzs3r1bfH19pVGjRtKsWTO5deuWqaU9kbSu7dmzZ2Xw4MFSvHhx6dWrlxw4cCBnhaVCdnwPYmNjpUWLFtK2bVsJDw9/Ytzo6Gjx8PCQ1atXJ4f5+fmJh4eHJCUlPRY/N3xv0e9Q0Dwr+fPnx8bGBgsLC7y9vfHz86NOnTo0adKE0NBQU8vLMHFxcXTu3JkqVarg7+/PyJEjcXZ2ZuHChbz33nvMmTOHK1euJMe/fPkyPXr04PXXX6d///7Mnj2b69evm/AMMsd///tfihYtytatWylRosQT444dO5bq1avz+uuvJ4ctXryY/v3757oBv3kdbW65GAsLC6ZNm4atrS1DhgzJcLpbt27xyy+/MGDAABo1aoSrqytFihTB1taWKlWq0KZNG/7zn/8k95tlNVZWVjRo0IB69eoRHR3NqVOnyJ8/P3Xq1KF69eocOnSIatWqMX36dBITExk+fDgODg507dqVF198kUOHDuHu7k779u1ZtWrVE5/Cmpo///yTVatWsXjx4nS7D/bv38/KlSuZN29eclhERAQbN26kT58+2S31+cPUVcdn2cytWfqAB1pDQ0PF29tbWrZs+cR3YN6/f1/8/Pxk7Nix4uXlJUWKFJH27dvLrFmzxM/PTwICAuTOnTty69YtOX78uPz2228yadIksbe3l1mzZkliYmKW6M0M586dE29vb6lXr564urqKj4/PQ8cjIyNl2bJl0qpVKylWrJi8/fbb4uvrK7du3ZJbt25JSEiIXLp0KdPaM6I1KipKgoKC5MyZM8lbQECAXLhwQcLDwyUmJkauXbsms2bNEkdHx+R3hz6JpKQkqVy5smzcuPGh8O+++066dOnyTHqzG/Jos9TkAp5lM0dzS0pKkm3btsl///tfsbOzk4kTJ0pCQkKqccPDw6Vz585StGhRqV27tnz88ceyZ88eiY2NzVBZZ8+elYYNG4q3t7eEhIQ8teanvbaJiYmycuVKeemllyQgICDNeJcvX5bp06dLrVq1xNbWVooVKybly5eXUqVKSZkyZWTmzJkSFRX1TFrPnj0rI0aMEE9PT7GxsZEKFSqIu7t78ubq6irlypWT4sWLS/78+cXOzk5effVVOXToUIbKDQoKkjJlyjzWr+bl5SVbtmzJtN6cJK+am35aakLu37/P+vXr+eOPP7h48SIXL17kwoULTJw4kbi4OP7991/Kly/PnTt3mDBhAuvWrUMpRdOmTenXrx/Ozs5s3LiRzp07U6lSJU6ePMmPP/7IlStX6N+/P4sWLXpi+W5ubvj5+fHVV19Rt25dpk6dSr9+/XKs78fCwoIePXqku+5cmTJlGDlyJCNHjnzs2LFjx5gwYQJjx46lcuXKuLu7Y2trS+HChZO3QoUKUbFixVSX+bl//z4jR45k5cqVvPPOOyxbtoxatWpl+RjD9evX07x584eurb+/P6GhoWY7LMbUZKm5KaVCgEggEUgQkbpKqQnAAOCGMdpoEfnNGH864A2MEJF9SqkKwHlgqIjMMcaZi+GXY2lWajU18fHxuLm54ezszKuvvoq3tzciwsWLFylatChnzpxh+PDhWFhY4OfnR+fOnfH19UVE2LNnD/369aN58+YcPnyYPXv2kJiYSP78+ZMntQ8ePDhDOiwtLRk1ahTt2rXjrbfeYtKkSXTu3JkOHTrg5uaGo6Njrn73ZY0aNdiwYQORkZGcPHmSoKAgoqKiiIqK4t69e1y/fp2oqChWrVpFz549+frrrwkLC6NixYqEh4fzySef4OrqSmBgIMWKFcsWjRs3bmTSpEn4+fk9FL548WLeeuutPDvXNtff71lZDQRCAPtHwiYAI1OJ6wFMBwoCa4xhFYBrQBCQ3xg2F3grtfLyerO0a9eu0rp1a+nWrZuULVtWSpYsKV27dpXVq1fL5s2bZe3atbJ69epUm2yRkZEybNgwKVWqlCxfvlwSExPlgw8+kFq1asmJEycy1Ex7tImUlJQkx44dk4kTJ0qDBg2kbNmyYmVlJcWKFZMaNWrIp59+KqdPn34sn9x4bVPj6tWr8ssvv0iXLl2kfv364uXlJRs2bEh1CEZWkJSUJJMnT5ayZcs+Ngzm2rVrYm9vL8HBwU/MIzdcW9Jolub0/Z7ZzZQ/GZZAEiBAynbQDeAPoC+Quxahz2K+++471q5dS8GCBZkyZQouLi4opdi7dy/Nmzd/YtrChQvzzTff0LNnTwYPHsyUKVPo1asXxYsXp3Pnzty8eZMmTZpQqFAhoqKiiIiIeGgLDw9n0KBBzJ07NzlPpVTymm0PplIlJSVx+/ZtAgIC+OWXX2jRogVubm4sXLgw1815TY+SJUtib2/P+vXrs72s2NhYBgwYwOnTp/nnn38oW7bsQ8fHjh3Lm2++iYuLS7ZrySXk+P2e1eYmgI9SSoDvROR7Y/hgpVQf4BCGKultETmplCoI/A48upTEVGCbUmpJFuvLVZQoUYJ33333mfKoX78+Bw4cwMfHhxUrVhAZGcn169fJnz8/QUFBVK5cmddee40KFSokDwcJDg7m9ddfp1GjRunmb2FhQYkSJWjYsCENGzZkxowZzJ8/n4YNG/Lpp59maojK88LFixfp3r075cuXx8/Pj4IFCz50/MiRI2zatIkzZ86YSGGWkavvd2WsCmZNZkqVEZEwpZQjsBMYAgQA4RguxOdAaRHpl0b6CsAWEammlFpmzKM+abTBq1evLhl9fV1UVFSembeXFVrj4+OJjY3l9u3bREZGUrVq1eRjcXFxnDlzhuLFi+Pg4JChSd6PEhcXx7lz5yhUqBB2dnbP1bVNi+joaK5evUpERASlSpWiVKlSj8VJSEjgzJkzlClTJkNzUHPD99bb2/uwiNR9NDyn7/dMkxVt2zTa4xN4pO2NoY3t/4Q0yccxtNH9gW8x0z63tMhqrc2bN5cVK1YIIOXLl5f9+/fL5cuXZcKECVK2bFlp1KiRrF27NtN9TxEREdKqVStZtGiRxMTEPJW206dPy4wZM9Ic7pLVZMf3ICYmRkaNGiUlS5aUuXPnSkRERKrxoqOj5cUXX5QxY8ZkOO/c8L0lA0NBcuJ+z+yWZTMUlFKFlFJFHuwDrQF/pVTpFNG6GE8gXUTkDHAKeCWrND6vvPnmm/z6669UqlSJixcv0qdPHz766CMGDRrE+fPnad26Na+99hp//vlnpvItUqQIW7ZsAaBXr15P9fq/0NBQRo4cSb169Zg5cyYhISGZzsOUHDx4kNq1axMcHMzx48cZNGhQqku8JyUlJfexff755zmuMyIigtOnTxMdHZ0l+eWF+z0rp1+VBH5XSh0DDgBbRWQ7ME0pdUIpdRzDY+DhmcjzC6BcFmp8LmnXrh07d+6ka9eufP755/j7+1OyZEleeOEFvv32W1q2bAnwVKvA5s+fHxcXF27fvs2XX36Z6fQtW7akRo0auLq6cvLkSerVq0fdunWZPHlyru6TiouLY/To0XTo0IHx48fzyy+/PHHJp1GjRnH9+nWWLFnyzOMIRYSkpKTkz0lJSYSFhfH333+zevVqpk+fzuDBg+nYsSM1atSgePHilC5dmg4dOlCyZElatWrFsmXLiI+PfxYZuf9+z4rqn6k23SzNOHXq1JGXX35Z5s2blxx25MgR6dSpk9jZ2QkgXl5eT5X3nj17JCgoSOzs7CQsLCzT6f39/cXe3l5Onjwp8fHx4uvrK4MHD5YyZcpIlSpVZOzYsRIaGvpU2lLT+qxcvXpVXnzxRenYseMTp8U9YO7cuVK5cmW5efNmpst6VO+ZM2fEy8tLihcvLk2aNBEXFxfJnz+/ODo6Sr169aRr167y4YcfyjfffCPr16+Xw4cPy40bN5K7HO7evSvr1q2TFi1aiJOTk3zzzTfpDhsij85QMLmAZ9meR3M7cuSIfPPNN9K/f38ZMGCAbNu2Te7fv59unuPGjZPOnTtLxYoV5fbt2w8du3jx4mNhT6P3o48+kn79+j1VHsuWLZMSJUrIzJkzk/vfEhMT5a+//pKhQ4eKnZ2dfPjhh8+8/NOzfg/Onz8vzs7OMn78+Az1UZ46dSpD49nS4oHepKQkmTdvntjb28u3334roaGhsnv3bgkICJDo6OinyvvAgQPStWtXcXBwkMGDB8uWLVtS7fvU5qbNLctIS+vevXvF0dFR3n//ffn222/lq6++Ei8vL2nRokW6ef79999StWpVGTJkiHTq1CndGyIhIUHi4+MzpffOnTtSsmRJOXz4cIbSPUpAQIA0b95c6tatK0eOHHnoWFhYmPTr10+qVKnyWC0uMTFRBgwYIHfu3Mmw1qelU6dO8tlnn2U4/muvvSZffvnlU5e3Z88eiYiIkLZt20q9evXkzJkzT51XWpw9e1amTJkiXl5e4u7uLjt27HjouDY3bW5ZRmpao6KixMXFJXmS9fnz52XkyJFSv359yZ8/f7qLJCYkJIiDg4MEBARI165dpWTJkjJp0qTHakIXLlyQ3r17S6FChcTS0lIqVaokPXr0kK+++kr8/PxS/WVPqfe7776TZs2aPfWo/6SkJFmyZIk4OjrKyJEjH2syTZ06VSpUqPDQTZ6UlCQYhh6k2/R7lu/Brl27xMHBQW7cuCFbt26VgwcPPjH+oUOHpHTp0nLv3r2nLnPPnj3Srl07efvttzNUQ38WHizaYG9vL/7+/snh2ty0uWUZj2qNj4+X7t27S69evZLDevToIUWKFBFfX1+JjIx8LI/IyEjp3r279OzZU3766Sfx9fUVd3d3mTt3roiInDhxQvr06SMFChSQ4sWLS+nSpaVChQpStGhR6dixowwcOFBq1KghVlZWUrRoUSlbtqy4u7uLh4eHrFq16qGlhlLqTUhIkBdeeEHWrVv3TNfg2rVr0qtXL6lQoYJs27btoWOLFy8We3t7mThxomzfvl1+++03ad68uQAyffr05HgrV66U+vXrPzQ041m+B8eOHZOGDRuKjY2NuLq6SoMGDZ4Yv23btg/1cT4NmzZtkqpVq2a7saVk6dKl4ubmltxVoc1Nm1uWsWfPHklMTJSYmBiJj4+X/v37S+vWrR9qSoaEhEiFChVk3Lhxqa5pNm7cOClWrJi8/vrr0qFDB2nWrJkMHjxYbty48VC8mJgYuXnzpoSGhkpgYKDExMRIlSpVpGDBgrJjxw75+++/ZerUqdKyZUspUaKEVKlSRYoUKSLVqlWTtWvXSmJiouzZs0ciIyPl7t27cu/ePfn+++8z1FTOCNu3b5eKFStKz549H+q8P378uAwfPlxatmwprVq1kuHDh8uOHTvE1dU1+XqsWbMmuUbXr18/OXPmzDN9DxYsWCC9evWSDRs2yOTJk0UpleaYtn379kmFChUkLi7uqcsTMZhbvXr1sm3+a1oMHDhQOnToILGxsdrcTLGZq7lt3rxZKleuLPny5RMrKytp0aJFqjdRQECAAOLn5/fYsZMnT8qkSZOkdu3a8sILL8ipU6cyXL6fn5+0bdtWbGxspGrVqtKzZ0+ZOnWq/PLLLzJ69GgpU6ZMsmm0adNG5s+fLwUKFJAiRYqIjY2NlCpVSsaOHftM1yAl9+7dk1GjRomDg4MsWrQozRs9KSlJ6tatK6+//rpMnjxZJk6cKBUqVJBZs2bJF198ISVKlJBff/31qTScOHFCHB0dZfLkydK6dWsZNmyYrFixItUflqSkJGncuLH8+OOPT1VWSnbv3i0NGjSQLl26PNXT1qclLi5OOnfu/OBJuja3nN7MzdzOnz8v77zzjsyaNUt+/fVXSUpKSl548siRI9K2bVuZOnVqcvxhw4YJIMuWLRN/f385fvy4/PXXX7Jr1y45evSoiBhutAULFkjFihUzXYuIjY2Vf//9V5YuXSrDhw+X8uXLS6tWreSrr76Sxo0bCyBVqlSRX3/99ZlrKBnh33//lXr16knTpk1TXZ1ExDBr4quvvpL//ve/Mnr0aFm6dGnyg5Hx48fL8uXLM13uzZs3pW7dujJr1qwMxf/tt9/E09MzS2ZdPFh8dPjw4eLk5JTj3+MLFy5oczPFZi7mFhgYKP369RM7OzsZM2aM7Nq1K/lYfHy87Ny5UxwcHOT111+Xbt26JR+7deuWLFiwQNq3by+VK1eWqlWripeXl3h7e0u5cuXkk08+SY7r5eUlu3fvfiadcXFxsnz5cnn//fflyy+/lKCgIBHJ2WubkJAgs2bNkhIlSsjEiRPTXXX4/v37EhwcLD4+PtK5c2eZMWOGAI81z9Pi1KlTUq1aNfnwww8z1DRMTEyUWrVqydq1azOUf3qkvLbbtm2T0qVLy5gxY3K0D06bmza3p2LBggVSokQJGT9+fPKTyz179khYWJh07dpVChcuLFWrVpVNmzaJr6+vNG3aNEP5hoeHi7Ozs6xYsUIiIiKkYcOGj63fn1WY4tpevHhROnbsKB4eHo81y8PDw2XAgAFSsWJFyZ8/vzg7O4u3t7f0799fVqxYIdOnT0/THH7//Xf573//K0OGDJGXXnpJ7OzsZP78+Rnu81qzZo3UqVMny/rIHr22V69elbZt20rdunVlw4YNqT5Mymq0uWlzyzSbNm2SihUrJteARAz9aKtWrZJSpUrJuHHjHhqqERgYKLa2tjJkyBBZvny5nD179ok30YEDB8TNzU1sbGySO4ezA1Nd26SkJFm3bp2ULVtWXnrpJZk2bZpMmjRJSpUqJcOGDZOzZ88+1lxOT+vevXvFw8NDAHnnnXcyNbg5Pj5eKleuLNu3b3+Ks0md1PQ+ePdEkyZNxNLSUtasWZNl5aVGXjW3vLm+sZmwZMkSBgwYwLlz55gzZw5bt24lOjqayZMns2HDBl588cWH4ru6uuLj48P+/fvZuHEjn3zyCRERERQrVgwbG5uHNnd3d1q0aIGfn1/ycXNDKcWrr75K69at2blzJ7t27SJ//vxs3ryZunUfW6EnQzRr1oyTJ0+yYcMGJkyYQK9evVi6dCkODg7ppl2+fDklS5akdevWT1V2Wty7d4+AgADs7e0pWbIk+fLlo127dgQGBhIZGanfwZAWpnbXZ9nyes2tatWqAkiDBg1k0qRJ8u+//0pSUlKmtN64cUPOnTsnp06dkn///Vf++usv8fX1lVmzZkmnTp2kWLFiUq1aNZk1a1aGRvA/Dbnx2qZFZrTGxcXJiBEjpFq1auleu9jYWHF2dpb9+/c/o8L/ExgYKIsXLxZbW1t54YUXkpd9t7CwkCJFikiNGjXk8uXLWVZeWqBrbprM4u+fodVgnoi9vT329vaPhXt7ezN06FASExP5888/mT17NmPGjMHDw4PWrVvTr18/KlWqlGa+UVFRbNu2jXv37mFra8vLL7/8VIta5mXy58/P9OnTCQ0N5bvvvmPUqFFpxl24cCFVq1alcePGWVb++PHjadu2LcHBwcn/46SkJJRS+u30GUC/cd7MsbS0pEmTJvzyyy+Eh4fz9ddfExcXx4svvkiPHj0eWjonJf/88w//+c9/2Lt3L99++y1lypTh7bff5ujRozl7AiZGKcVHH33EvHnzSEhISDXOvXv3+OKLL5g0aVKWlr1v3z7s7Owe+vGysLDQxpZBtLk9R1hbW9OkSRNmzJiR/I7UJUseX7b++vXrXLx4EaUUP/zwA7t27eLkyZNUrVqV1q1b880332BorTwf1KlTBycnJzZu3Jjq8Tlz5tC0aVNq1aqVpeWOHTuWCxcucPjw4SzN93lBN0ufUwoUKMD8+fNp3bo158+fJzExkYsXL3LgwAHCw8OpX78+n3zySXIt4cGLkbt27Ur37t25efOmSVaUNRXDhg1j1qxZvPbaaw+F37lzh6+++or9+/dneZnvvvsua9eupWPHjri5uVGtWjUqVaqEq6tr8mZlZZXl5ZoLuub2HFOzZk327t1LbGwsRYsWpW3btvz666/cunWLHTt28NFHj76kCCpWrMjKlStZsGABsbGxJlBtGjp37kxISAgHDhx4KHzUqFF06tQJDw+PLC9TKYWDgwNBQUF89NFHuLm5cf78eebPn0+nTp2wt7enU6dOzJ07l7Nnzz5XtemMoGtuzzlVqlThq6++ylSaSpUqUatWLdavX88bb7yRTcpyF1ZWVkydOpWOHTvSo0cPEhMTOXr0KLdu3eKvv/7K1rILFChA+/btHwu/ceMGu3btYufOnUydOhULCws8PT2pWLFi8tawYUPKlXs+V+rX5qZ5Kt577z0mTJhAly5dTC0lx+jVqxcuLi78+eefWFlZ0bZtW1q2bGmyMYQODg707NmTnj17IiIEBQURGBjI+fPnOX/+PP/88w/vvfcederUoW/fvnTp0oVChQqZRKsp0OameSq6dOnCli1b8Pb2ZuLEiaaWk2M0aNCABg0amFrGYyilcHNzw83N7aHwmJgYfv31V3788UeGDBlCly5d6Nu3L02aNMHCwrx7pcz77DTZgo+PDx9//DGdO3fmxRdf5NSpUyxbtszUsjSpUKBAAXr06MG2bds4deoUVapUYfDgwVSqVIlPP/2UoKAgU0vMNrS5aTLE1atXWb16Ne3bt+eNN97gwIEDzJo1iw0bNlCwYEEqVqxoaomadChdujQjR47k+PHjrFu3joiICBo2bEiTJk3Yvn272T2Q0M1STbqICC4uLtjZ2VG3bl0GDx7MrVu3iI2NpXHjxri7u1OzZk1u3bqFnZ2dqeVq0kEpRe3atalduzbTp09n/fr1fPjhhxQrVozPPvuMli1bmsVAYV1z06SLiFC/fn1sbW2JiYkhMjKSSpUqUbduXU6ePMnJkycpU6YM5cuXx9fX19RyNZnAysqK7t27c+LECYYMGcKgQYNo3rw5+/btM7W0Z0bX3DTpYmFhwZ49e1I9NmDAALZs2UJISAiNGzcmMTExh9VpsgJLS0t69uxJt27dWLlyJe+88w42Njb06NHD1NKeGl1z0zwTSimKFClCwYIFuXTpUpZPQdLkLPny5ePNN98kICCA+fPnc+XKFVNLemq0uWmyhAIFCphNc0ZjqK03btyYuXPnmlrKU6PNTZNlFC9enOjoaFPL0GiAPG5uMTExppagSUFUVNRzt+abJveSp80tMDCQH374wdQyNEaCg4OfuACmRpOT5Glzq1SpEmPGjOHEiROmlvLcIyLa3DS5ijw9FKRQoUKMHz8eb29vRo8ezdChQ8mXL0+fUq5ADcjcAM4ZjWfg3csbkqDoiKLZpCpjyELzGmWveXrydM0NoEePHvz5559s2bKF+vXr61VLTcVloJSpRWg0/yfPmxuAu7s7u3fvZujQobRr147hw4cTFRVlalnPF2eBCqYWodH8nwy14ZRSIUAkkAgkiEjdFMdGAtMBBxEJz2hapVQZYLnxWC8RiVJKTQBGARVE5LoxXpSIFM6ARvr27Uv79u0ZOXIkVatWZe7cuXTo0CEjp6hJg4w089avX49bETdO7Txlkm6BzDajzYmoqCiuXLlC/vz5cXZ2ztGyc7svZOab6P2oSKWUE9AKuJjZtMBQYAjgAvQGFhjDw4ERwH/TExQSEkL37t2JjY0lPj4eJycnPDw86NatG23btuXDDz9k+fLlqb4ERZM1LF++nJs3b7J+/Xrd35kDhIWFMX78ePz8/AgLCyMxMZHSpUtz7949bGxsaN68Oc2bN8fb2zunzC7X+cIDnrVZOhODoz5NL64lkGTcUv70LgG6K6XSXV6icOHCdO7cmbfffpuBAwdSvXp1zp07xyuvvML06dM5ceIEhQoVokWLFmm+lk3zdJw4cYL27dszceJEKleuTLVq1UwtyayJiYlh0qRJVK9enRIlSrBx40ZCQ0O5d+8ewcHBXLlyhe3bt1O/fn22bduGl5cXFStWpF+/fixfvjyn33dhUl94QEbNTQAfpdRhpdR/AJRSHYHLInIss2mNzAW+A97DUA19QBSGExmWnih7e3t69uxJ586deeWVVxg0aBBvvfUWxYsXZ8KECdjY2LB48WJatWpFQEAA9+7dy+DpalLj/v377Nixg969e/PSSy/RunVr/P39TbbM9vNCeHg4L7zwAseOHePAgQNMnToVT09PihYtmrw0kVIKDw8P3n//fVavXs3Vq1fZunUrdevWZfny5bi6ujJv3jzu37+fldJypS88QGVkgTqlVBkRCVNKOQI7MVQbpwOtReSusf1cN4229WNpRcQvjXImGE9iEXAUqA6EpdW2rl69usyePTv58/Xr17l69SplypR57C3sN27cwNbWFmtr63TP19RERUVRuHC63YzZxuEL/3/iXMWxCtHR0URGRnL37l2sra0pXrw49vb2WFpaAqnrFRHi4uKSt9jYWBISErCxsaFYsWIULFgwW/TWca7zxLimvraZJSoqimvXrmFtbf1ML3qJjo7m8uXLJCUl4ebmlqklxr29vQ+n7E97QG71hQdkqJNERMKMf68rpTYAzYCKwDHjL0c54IhSyktErqaT1gtI9SRSpLmjlPoZGPikePnz56d58+YkJCQwduxY1q5dy+7du1Pta5g7dy61a9fOE4NM9+7dS/PmzXO0zJUrV/LHH39w584dVvitgGjgDni4eVCzZk0aNWpE165dH7rBIiMjCQ4O5vbt2/zzzz8EBQURFBREcHAw169fx9nZ+aH3bDo4OHDixAlWrFhBtWrVmDhxInXrPnbPZBrvAd7J+9L3yT/Wpri2z8K6deuYNWsWBw8efOYf5qSkJIYMGcLKlSvZtm3bM0+Vy62+8IB0zU0pVQiwEJFI435r4DMRcUwRJ4RUHDqttBkRBnwNHHySxjt37jBq1Cg2btyIs7Mzf//992M1NoBz586RmJhI2bJlM1j088ecOXP466+/aNSoETgDdkBx+OPrP5IN64cffnjIwCIiIqhUqRIffPABN27coE6dOnTv3p1KlSrh5OSU6gOGnj17Mn78eBYvXkznzp2pW7cuEydOpEaNGjl+zrmdixcvEhoayk8//ZQlLQ4LCwvmzJlD79696dmzJ2vXrn3qh0C52RcekJG6aUngd6XUMeAAsFVEtqcVWSlVRin129OkTYnxgmwA0vyvhoeHY2Njw4oVK9ixY8dDxiYi+Pv7M2nSJJo2bYqjo6PuG3oC27dvZ/Xq1YZa70FgH7DK8BLm999/n40bNxIbG0uzZs2YMmUKR44c4d69e5w4cYJKlSoxY8YM3nvvPV566SUqVqz4xJvG2tqagQMHEhgYiLe3N02aNGHkyJE5dq55hT/++IP8+fOzefNmLly4QHx8PAEBAfj7+xMQEPBUYzktLCxYunQp0dHRDBw48Fnem5BrfSFl5Dy71alTRx7l6NGj8u6774qTk5M4OzvL+++/L//++6/s2bPnsbi5FVNrpR9CF4TeSFJSUrrxn0Vv//79pVKlSuLj45Nu3NjYWDl37pxcu3ZNoqOjk7XxDslbdmrNaUJCQmT+/PnSqVMnqVGjhpQrV05cXFykSpUq4ubmJgULFpQBAwY8Vd4RERFSp04d+fTTT9ONCxySXHC/Z3Yzm4FJp0+fZsKECfj5+TF06FB27NiBh4dH8tOkvXv3mlZgXsICKGHYze4XhWzfvh0fHx+qVKmSZpyEhAQ++ugjvv/+e0qUKJH8HoeEhAQqV65sGDzg+n/N5oKzszMeHh7079+f2bNnU716dVq1apV8/N69exQrVoxZs2Zluv+sSJEi/PbbbzRq1IjSpUvz3nvvZbV8k5Pnzc3Pz4/vv/8eHx8fRowYwZIlS56rt2rndXr37s3kyZNZvnz5Q+Hx8fHcvXuXM2fOMGXKFBITEwkJCcHBwSE5zv379zl69Cj136sPOwA7w/ehSZMmZvH2pgdYWVkxYsSIx8JFBEtLy1T7444fP87q1au5du0apUuXTt7KlClD6dKlKVWqFI6Ojmzfvp0mTZrg6OjIq6++mhOnk2PkaXM7duwYgwYNok+fPsybN4+iRU27IoUm84wePZpy5cpx9+5dgoOD+eCDDzhy5AixsbEULVqU8uXL0717dz788EPy58//UNr8+fPj5eUFdYAaQBD0798fBwcHPv74Y1555RWzfqu6jY0NL774IvXr16dWrVqUKlWKpKQk7t69y9y5cylWrBhffvklV65c4fjx4+zYsYMrV64QFhbG9evXsbW1pXTp0hQoUCD5XbTVq1c39WllGXna3Dw9PTl27JhZ/Uo/bxQoUIACBQpw/vx5OnTowOeff87WrVspXLhw5v6v+QAPOLP3DOvXr2fixIkMHTqUl156iRYtWtCiRQtKlTKvZUvy5cvHrl272L59O5cuXeLatWtYWVlRoUIFPvzwQ958801q1qyZatqkpCTCw8O5cuVK8layZMmcPYFsJk+bW/78+bWx5WHi4+N5//33qVmzJkeOHKFOnTr069fvmfK0tLSkW7duvPbaa5w6dQpfX1/WrFnDoEGDKFOmDB9//DEbN26kefPmFCtWLGtOJJu4f/8+8fHx3Llz5zGtd+/eZf369WzcuJGyZcvSoEEDOnbsmOHhThYWFjg6OuLo6Gi2w3DMt86uybWEh4fzyy+/UK9ePcLCwli7di3ff/897777bpaVoZSiatWqDBkyhI0bNxIeHs7SpUuxsrLi22+/xcnJCS8vLz7++GN8fHxM8mKbmJgYfH19GT9+PL169aJt27bUrVsXpRRKKaytrTl+/DhNmjQBDGa3efNmunfvTvny5dm0aROvv/46RYoUoU+fPgwaNCjHzyE3o81Nk6MsXLgQBwcH5s+fn9wE/euvv7h69Spt27bNtnItLS2pV68epUqVwsfHh/DwcKZPn461tTWfffYZjo6ONGvWjM8++4xz585lm47Y2FhWrVpF27ZtcXBwYMyYMcTFxdG2bVuGDRvGt99+i7u7+0Np/P396dChA2XLlmXatGm0aNGC8+fPs379ehITE/nhhx8YPXo0q1evzjbdeZE83SzV5D3i4uKwsbHht99+w8bGhrCwMPr27cvPP/+cPFc1J7C2tqZZs2Y0a9aMiRMnEhUVxf79+/Hx8cHLy4s+ffowduxY7OwyvAjFE7l79y7Tpk1jwYIF1K5dm7fffpu1a9emOs81ICAgeX/v3r1cv36dbdu2MWvWLJycnHB1deX8+fMEBgYSGBjI9u3bqV27dpboNCd0zU2TM6xYARUqMGjoUM6LsKNvX+Lj4+nRoweDBg3C29s7/TyykcKFC/Pyyy8zc+ZMTp48SUxMDJUrV2bGjBnPtFzQ/fv3mT17Nu7u7oSFhXH48GF27tzJG2+8keEJ/A4ODvTp0wcXFxfy5ctHdHQ0X375Jba2tvz111/a2NJAm5sm+1mxAv7zH7hwASVCqbg42qxbx2ceHjg4ODB69GhTK3yIkiVLMn/+fPbv38/+/fvx9PRk5cqVJCUlZTgPEWHNmjV4enqybds2du7cyQ8//ECFChWeSZtSir///hsHBwc9pjMdzM7cEhISuHTpEhcvprcIqCbHGDMGHumwt0lM5KPbt1m7dm2uHYvm4eHBr7/+ytKlS5k5cyb169dnz549TzS5wMBAxo8fj6urK19++SXfffcd27Zty9LxY5UqVcLGxobAwMAsy9McMZs+tytXrjBlyhQWL15MfHw8NWvW5MCBA6aWpQFI44fG9s4dyANDeZo1a8bff//N6tWref/997l48SIuLi64u7tjZ2dHWFgYoaGhXL58mXz58tGzZ0/WrFlD7dq1s22oUqNGjfjjjz8ee/ig+T952twOHz7MvHnzCA4OZunSpbz99tts376d119/na+//trU8jQPKF8eLlxIPTyPYGFhQc+ePenZsyf37t0jKCiIs2fPcvv2bcqUKUO5cuUoW7Ys9vb2OTL2skmTJvj6+vL2229ne1l5lTxtbgBffPEFPXr04OTJk9y5c4fWrVvzxRdf0LhxY1NL0zzgiy8MfW4pm6YFCxrC8yCFChWiRo0aJh382q1bN8aMGcPNmzcpUcLMVgzIInJnZ0cGqVGjBmFhYcyYMQNfX1+aNWvG5MmTn3mUuyaL6dULvv8enJ0NzVBnZ8PnXr1MrSzPUqhQIWxsbLh8+bKppeRa8nTNLSgoCEdHR27fvk3t2rXZsmWLYSK1JvfRq5c2syzEx8cHT09Ps5rontXkaXMrW7Yse/fuxc7OLk+8+EWjySry5cunv/PpkKfNrUiRIpQuXdrUMjSaHMfT05MTJ04gInrxiDTI031uGs3ziouLC4mJiXo85xPQ5qbR5EGUUjg7O3P+/HlTS8m1aHPTaPIgly5dwt/fP0ve+2quaHPTaPIY+/bto1GjRnz++ecZnnz/PJKnHyhoNM8TiYmJ+Pj40LdvX3744Qfat29vakm5Gm1uGk0uJygoiB9++IEff/yRsmXL8s0332hjywDa3DSaXMi9e/dYu3YtsbGxdO/end69e7N9+3aqVatmaml5BrMwt9u3b/P7778TEhJCaGgo4eHhVKpUiVq1atGmTZtcu6SORpMSEeHvv/9myZIlrF27lsaNG9O/f38uXbr02GsNNemTp+/6a9eu0bhxY5ydnZk7dy5nzpyhWLFi1K9fn7t379KuXTu95pUmT7Bz506qVq3KW2+9haurKydPnmTz5s0UK1ZMG9tTkqdrbhEREYwbN46mTZtSoECBh47Fx8fz9ddf69HbmlxNYmIiEydOZPHixSxatIi2bdvq72wWkafNzc3NjTZt2qR6zMrKiqlTp9KkSRNGjx6tJxhrch1Xr17ljTfeQCnF4cOHze6l0aYmTzdL02PEiBHs3LmTXbt2cerUKT2aW5NjiAjh4eEkJCQkhyUmJnL69Gl++ukn+vTpg4eHB02aNMHHx0cbWzZg1uYGUL16dTZv3oydnR1Dhw41tRzNc0BwcDD169fHxcUFNzc3+vfvT61atShevDivvPIKW7ZsoV69egQGBjJx4sQcfaXh80SebpZmhsKFC3Pv3j1Ty9CYOceOHaNdu3Z88sknDBw4ED8/P/744w/effddXF1ds+w9qJr0eW7MLTo6Wo8RegrUgPQ7t2c0noH3ANO+dzQ3sH//fl577TXmzp1Lt27dAGjevDnNmzc3rbDnFLNvlj7g7t27tGjRwtQyNGbK6tWr6dq1KytWrEg2No1pydM1t8DAQDp16kSVKlXw8vKiffv2qY4JCgwMJCYmhpdfftkEKjXmzpo1axg+fDjbt2/Xb3/PRWSpuSmlQoBIIBFIEJG6KY6NBKYDDiISbgybDngDI0Rkn1KqAnAeGCoic4xx5gKHRGTpo+U5OjrSp08fTp8+zaxZs3jvvfd46623eOedd3BzcwMgNjaWTz75hFdffVUvy5xBZKFkKv7evXuRvplLY06UKVMGS0tLXFxcTC0lR8np+z2zZEez1FtEaj5yok5AK+BiijAP425TYFCK9NeBYUqpdIdlFy1alK5duzJ27Fj27t3L/v37EREaN25MuXLlaNq0KY6Ojty7d08vR67JNho3bkynTp0YMGAAiYmJJtUiIkRERHD+/HliYmJyosgcu98zS041S2cCo4BfU4RZAkmAACl7rW8AfwB9gYWZKcTd3Z1p06bx5Zdfcv78eYKCgqhbty52dnbs3bv3mU5Ao3kS06dPp0OHDvTu3Zt58+Zly1PRCxcu8Oeff3Ljxg3Cw8PT3KytrSlevDh37txh4cKFdO/ePcu1pEOO3O/pkdXmJoCPUkqA70Tke6VUR+CyiBxLOa1ERE4qpQoCvwMfPZLPVGCbUmrJ04hQSuHi4vLcNRM0pqNAgQJs3ryZDz74ADc3N1577TXq1q1LlSpViI+P586dO9y5c4dy5crRokWLDC3mICLExsYyadIkNmzYwMWLF2nevDmlSpXC3t6eqlWrYm9v/9BWokQJbGxsADh8+DAvv/wyr7zyCoUKFXqq80rnvai54n5PCyWSdX0lSqkyIhKmlHIEdgJDMLS7W4vIXWMbve6DNngq6SsAW0SkmlJqmTGP+qTRBq9atarMmzcvQ9qioqLyzKqleUkr5C29OaH1/v373Llzh+joaGJjY7GwsMDS0hILCwuio6Oxs7NLs5vk7t27REREEBsbS2xsLCVLliQuLo7ixYs/le6QkBDi4uIoWrQoNjY2FCxYMFMT8YOCghgwYMDhlM3OB+T0/Z5pRCRbNmACMA5DmzrEuCVgaIeXSiNNBcDfuO8B+APfAm+lFr9w4cISFxcnGWHPnj0ZipcbyEtaRfKWXlNrDQgIkNKlSz8WfuzYMXnxxRelXr16MnXqVPntt9/k9OnTz6z3/v37smzZMvnoo4/klVdeETs7O9m3b1+G0sbGxgqG2tkhyQX3e2a3LGuWKqUKARYiEmncbw18JiKOKeKE8AQnT4mInFFKnQJeAQ6kFsfS0pLBgwfz3Xff6ZUUNHmCihUrEhsby08//YS7uzsiwrJly1i7di1ffPEF/fv3f6jJevXq1Wcqz8rKijfffDP5c+vWrYmIiEgz/rlz5yhXrhz58+fH2tqao0ePUrNmzcfimeJ+zyxZ2edWEthgNJl8wM8isv0Z8/wC+Detgy4uLvz999/Mnj2bYcOGPWNRGk32Y2Vlxffff8+KFSu4fPkyCQkJtGnTBn9/fxwdHdPP4Bn5999/qVWr1mPhhw4dYvz48fz111+UKlWKf/75hyJFilCjRo20ssrx+z2zZJm5icg5IM0rYYxTIZ3jIUC1FJ+P8YThKhYWFmzatIkGDRrg7u6uB+lq8gSvvfYar732mknKLlOmDOPGjaN27doULFgQCwsLNmzYwOHDhxk9ejTr169n8ODBdO3albZt26aZjynu98yS56dfVahQgbVr19KnTx++/fbbB215jUaTClu2bMHDw4PTp0+zf/9+fH19adOmDUFBQQwcOBBra2vmzJlDmzZtCA0NJTQ01NSSn5o8Pf3qAY0aNeKPP/6ge/fu7N69m0WLFlG8eHFTy9Joch1OTk6MGjXqiXFsbGwYMWJE8ueZM2dmt6xsIc/X3B7g7u7O33//Tbly5fDw8GDixIncuHHjoTi7du1i2rRpbNy4keDgYF3L02jMGLOouT3A2to6eY7pzJkzcXd3p0OHDjRo0IACBQowbdo0WrdujZ+fH//++y9WVla0bt2aVq1aUbZsWSwtLbG0tCRfvnxp/s2XLx8lSpQgXz6zunQajdlhlneop6cn33//PZ9//jnr1q3j8OHDNG7cmAMHDiQPhBQRzpw5g4+PD8uWLePmzZskJCSQmJhIYmJiqvsJCQkkJCRw584dnJyccHV1xcvLi969e+Pu7m7is9ZoNCkxS3N7QMmSJRk4cCBgWLki5QhvpRSenp54enpmaBhJYmIi27dvp0CBAjRq1IgLFy4QGBjI7t27adasGeXLl+ftt9/mrbfeSp7+otFoTIdZm1taJCYmcunSJQIDAzl79iyBgYEEBgZy/vx54uPjUUphYWGBUip5/9atW5QrV4579+5x8+ZNypUrR1JSEnFxcRQrVoyCBQuyadMmvvjiC8aMGUO/fv30+yY1GhPyXJmbv78/7733HocPH8be3h43Nzfc3Nxwd3enRYsWuLi4YG1tjYiQlJT00F8bGxtcXV0REUJDQ7l69SoWFhZYW1sTHR3NJ598Qnx8PAsWLGD27Nl8+eWXjBs3jj59+uj+OY3GBJjFXRcfH09ISAhnz55NromdPXuWkJAQrKysKFy4MAMGDKBbt25MmTIFHx8fChYs+FRlKaVwcnLCycnpoXAfHx/Gjh3LBx98gJ+fH8HBwYwbN44pU6bw5ptv0rhxY+rXr//UqzNoNJrMkafN7ejRo9ja2hIXF0fZsmVxd3fHzc0NT09POnXqRMWKFUlMTCQqKoqrV6/y77//Uq5cuWzRYmlpyZQpUyhcuDCtWrVi79697NmzBz8/P7Zu3cq4ceM4evQo8+fPp0+fPtmiQaPR/J88bW7VqlVj3759FChQIN3+rb1792absaVk9OjRREVF0aZNG3x9fWnatClNmzYFYPbs2Rw8eFCbm0aTA+TpQbz58uWjaNGiuarjXinF5MmTadiwIe3atXvoXakVK1bk4MGDevCwRpMD5Glzy60opZg1axaVK1emU6dOxMbGAtCuXTtu377Nnj17TKxQozF/tLllExYWFixcuBB7e3u6detGfHw8lpaWjB49ms8++0zX3jSabEabWzZiaWnJTz/9hFKK3r17k5iYSK9evbh9+zYLFiwwtTyNxqzR5pbNWFlZsWbNGm7evMk777yDhYUFa9euZfz48Rw4kCULjmo0mlTQ5pYD2NjY8Ouvv3L27FmGDRuGq6srCxcupGvXrmzfvl03UTWabECbWw5RqFAhfvvtNw4dOsQrr7xC/fr1mTdvHsOGDaNly5bMnz+fwMBAU8vUaMwGbW45SNGiRfHz86N27drUrFmTpKQk/P39efvttzlw4ABNmzalRo0aXLlyhYCAAFPL1WjyNNrcchgrKys+//xz1q9fz4gRI+jbty/16tXjhx9+IDQ0lLlz55KQkECLFi1o3LgxQUFBppas0eRJtLmZiIYNG3L06FGqVKlC06ZN6dmzJzdu3KBJkyY4OTlx8eJFatasqZ+qajRPiTY3E1KkSBHGjh1LcHAwFStWpE6dOpw9exYwDCPp3LkzBw8eNLFKjSZvos0tF1CkSBEmT57Me++9xxdffJEc7uLiwsWLF02oTKPJu2hzy0UMGjSIzZs3Ex8fDxjeMRkWFkZiYqKJlWk0eQ9tbrkIOzs73n77bUJCQoiKisLGxoYXXngBX19fU0vTaPIc2txyGVOnTiV//vw0b96cq1ev0qdPH5YtW2ZqWRpNnkObWy7DysoKZ2dnOnbsSMOGDalVqxabN29m586dppam0eQp8vRilebMp59+irOzM6+99hpDhw6lf//+tGvXjunTp1OkSBFTy9Nocj265paL6du3L8uXL2fBggWMHTuW+Ph4PD09mThxon6KqtGkgza3XE6rVq3YuXMnn3/+OcWKFWPZsmXcuHGDWrVq0aZNG9asWUNcXJypZWo0uQ5tbnmAGjVqcPDgQaKioujUqROnT5+ma9euhISE0L17d4oWLcqwYcPYvn07ERERppar0eQKdJ9bHqFUqVJ89913TJ8+PfnVgaVKleLgwYNs376d2bNns2zZMuLj43F3d8fV1RU7Ozvs7OwoXrz4Q/v29vZ4enpiaWlp6tPSaLINbW55DFtbW1555ZWHwq5cucKiRYv4/vvvcXR0xMvLi/r16xMTE8OtW7e4du0ap0+f5tatW9y+fZsrV65w7949unfvzhtvvEGdOnVQSpnojDSa7EGbmxlQunRpxo0bxyeffMLWrVuZP38+S5YsoVy5cri4uODi4kLlypWT911cXLh69SorV66ke/fuNG7cmB9++AELC91LoTEfMmRuSqkQIBJIBBJEpK5S6nOgE5AEXAfeEpGwjKQ1hpcBlhuP9RKRKKXUBGAUUEFErhvjRYlI4Wc5yeeFfPny0alTJzp16kRcXBwXLlzg3Llzyds///xDcHAw586dw8rKChcXFzw8PFi2bBmtWrWid+/epj4FTR4it/tCZmpu3iISnuLzdBEZZyxoKPAp8F4G0wIMBYYALkBv4MHaPuHACOC/mdCmeQRra2vc3d1xd3d/7JiIcPPmzWTTGzVqFA0bNjSBSo0ZkGt94ambpSKS8rFcISCzLwKwxODuSUDKDp8lwFtKqS9F5NbT6tOkjVIKe3t77O3t8fLyMrUcjRmRm3who50sAvgopQ4rpf7zIFAp9YVS6hLQC4NDZzgtMBf4DoOrL08RHmU8kWEZ1KbRaExDrvYFlZE3LymlyohImFLKEdgJDBERvxTHPwFsRGR8ZtM+EneC8SQWAUeB6kBYWm3r6tWry+zZs9PVDxAVFUXhwnmj6y4vaYW8pTcvaYXcodfb2/vwgz6xlORWX3hAhpqlDzoEReS6UmoD4AWkFPIzsBV47CQykDa18u4opX4GBj4p3oPVMzLC3r17MxzX1OQlrZC39OYlrZC79eZWX3hAus1SpVQhpVSRB/tAa8BfKeWWIlpH4ExG02ZEGPA18C56uIpGk+vIC76QEeMoCWwwDvLMB/wsItuVUuuUUpUxdPxdwPhExPgod5GItEsrbUbOQETCjY4+PCPxNRpNjpLrfSFdcxORc0CNVMK7phE/DGj3pLRPKGvCI58/BD7MaHqNRpMz5AVf0EPSNRqNWaLNTaPRmCXa3DQajVmizU2j0Zgl2tw0Go1Zos1No9GYJdrcNBqNWaLNTaPRmCXa3DQajVmizU2j0Zgl2tw0Go1Zos1No9GYJdrcNBqNWaLNTaPRmCXa3DQajVmizU2j0Zgl2tw0Go1Zos1No9GYJdrcNBqNWaLNTaPRmCXa3DQajVmizU2j0Zgl2tw0Go1Zos1No9GYJdrcNBqNWaLNTaPRmCXa3DQajVmizU2j0Zgl2tw0Go1Zos1No9GYJdrcNBqNWaLNTaPRmCXa3DQajVmiRMTUGp4apdQhU2vQaJ4DwkWkralFZJY8bW4ajUaTFrpZqtFozBJtbhqNxizR5qbRaMwSbW4ajcYsMRtzU0oNU0r5K6VOKqU+MIaVUUr5KqV+VUoVNoZNUEpdVkodTbEVy2ltjxxvrpS6m0LPpymO9VBKHUmZTikVopQ6kSL+7GzQvEQpdV0p5Z8izE4ptVMpFWj8WzzFselKqUNKqWbGzxWUUjGPXOc+Oa3rkbQpr9uhFOE5+j1J4xy6Gb8fSUqpuo/Ez5Fra3aISJ7fgGqAP1AQyAfsAtyAqUBVoAPwnjHuBGCkqbU9Eqc5sCWN9BsBS2AVUNgYFgLYZ7PupkBtwD9F2DTgY+P+x8CXxn0PYLrxHNcYwyqkTGsKXamkTfW65fT3JI1z8AQqA3uBuinCc+zamttmLjU3T+BvEYkWkQRgH9AFgykkGTeVy7RllAe6hRw8BxHxA249EtwJ+NG4/yPQ2bj/4Dpnu8ZM6sooOfo9Se0cROS0iAQ8QVuO/v/NAXMxN3+gqVKqhFKqINAOcALmAt8B7wHLU8QfnqI6v8dE2h6lgVLqmFJqm1Kqaorw9cAh4JCIRKYI35PiHIZnn/yHKCkiVwCMfx2N+ycx1Cx+B+aniF/pkaZTk5zUlQoC+CilDiul/pMiPDd8T1IlF1zbPEs+UwvICkTktFLqS2AnEAUcAxJE5AKGJsCjzBSRGabU9ki0I4CziEQppdphaIq6GdP/yP9rJSnxFpHwbBOeSURkSCrBwSJSM6e1PIFGIhKmlHIEdiqlzoiIX274njyJPHJtcx3mUnNDRBaLSG0RaYqhyh9oak0PSE+biESISJRx/zfASillbwKp6XFNKVUawPj3uon1PCBDukQkzPj3OrAB8MoxhZocx2zMzfhrjFKqPPAqsNK0iv5PetqUUqWUUsq474Xh/3Izp3VmgE1AX+N+X+BXE2pJSbq6lFKFlFJFHuwDrTF0GWjMFVM/0ciqDdgPnMLQ7Gv5hHgTgMvA0RRbhZzWhqF/58GTucHASePxv4GG6eQXApxIoX9ZNmheCVwB4oFQoD9QAtiNoea5G7B7QvoKQMwj13loTuoCygC/GfddjNf3mPFaj0mnnGz7nqRxDl2M+3HANWBHTl9bc9v0xHmNRmOWmE2zVKPRaFKizU2j0Zgl2tw0Go1Zos1No9GYJdrcNBqNWaLNTaPRmCXa3DQajVnyPwweCdm1c7uFAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import matplotlib.pylab as plt\n", "from matplotlib.patches import Rectangle\n", "import cartopy.crs as ccrs\n", "\n", "ax = plt.axes(projection=ccrs.PlateCarree(central_longitude=0))\n", "ax.coastlines()\n", "ax.set_extent([left_bound - 1, right_bound + 1, bottom_bound - 1, top_bound + 1])\n", "ax.plot(max_plancks_birthplace_x, max_plancks_birthplace_y, \"-ro\")\n", "ax.add_patch(\n", " Rectangle(\n", " (left_bound, bottom_bound),\n", " np.abs(np.diff([left_bound, right_bound])),\n", " np.abs(np.diff([bottom_bound, top_bound])),\n", " edgecolor=\"darkgreen\",\n", " facecolor=\"none\",\n", " lw=3,\n", " )\n", ")\n", "gl = ax.gridlines(draw_labels=True, x_inline=False, y_inline=False)\n", "ax.set_title(\"Max Planck's birthplace\")\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "... that looks good ! Now we select the corresponding cells within the dark green window and display the indices." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Cells" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([4282376, 4282377, 4282378, 4282400, 4282401, 4282402, 4282403,\n", " 4282404, 4282405, 4282407, 4282408, 4282409, 4282410, 4282411,\n", " 4282412, 4282413, 4282414, 4282415, 4282420, 4282421, 4282422,\n", " 4282480, 4282481, 4282482, 4282483, 4282484, 4282485, 4282486,\n", " 4282487, 4282488, 4282489, 4282490, 4282491, 4282493, 4282497,\n", " 4282501, 4282508, 4282510, 4282511, 4282512, 4282513, 4282514,\n", " 4282515, 4282516, 4282517, 4282518, 4282519, 4282520, 4282521,\n", " 4282522, 4282523, 4282527, 4282558, 4282936, 4282938, 4282939,\n", " 4282941, 4282960, 4282961, 4282962, 4282967, 4282968, 4282969,\n", " 4282970, 4282971, 4282972, 4282973, 4283128, 4283129, 4283130])" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "window_cell = (\n", " (grid.clat >= np.deg2rad(bottom_bound))\n", " & (grid.clat <= np.deg2rad(top_bound))\n", " & (grid.clon >= np.deg2rad(left_bound))\n", " & (grid.clon <= np.deg2rad(right_bound))\n", ").values\n", "\n", "(window_cell_indices,) = np.where(window_cell)\n", "window_cell_indices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "What we have shown above are all the indices of the cells within the green window of the native grid. Since we selected them via `np.where()` we obtained the indices in 0-based python thinking. We now do the same for the vertices and edges that appear in our selected window:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Vertices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We select with the function `.isel()` and the `vertex_of_cell` information from the grid all the vertices of the cells we have cut out one step before. We also sort the indices of the vertices and delete all duplicates via`np.unique()`. Since the ICON code is written in Fortran, i.e. the integer-values are 1-based, we subtract `1` to get back to python thinking, hence 0-based." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([2144929, 2144932, 2144933, 2144935, 2144937, 2144938, 2144939,\n", " 2144953, 2144954, 2144955, 2144956, 2144957, 2144958, 2144959,\n", " 2144960, 2144961, 2144962, 2144963, 2144968, 2144978, 2144983,\n", " 2144984, 2145002, 2145003, 2145004, 2145005, 2145006, 2145007,\n", " 2145008, 2145010, 2145011, 2145014, 2145020, 2145021, 2145022,\n", " 2145023, 2145024, 2145025, 2145026, 2145217, 2145218, 2145220,\n", " 2145222, 2145223, 2145224, 2145231, 2145237, 2145238, 2145239,\n", " 2145286], dtype=int32)" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "window_vertex_indices = (\n", " np.unique(grid.vertex_of_cell.isel(cell=window_cell_indices).values) - 1\n", ")\n", "window_vertex_indices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Edges" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Same as for the vertices, we select with the function `.isel()` and the `edge_of_cell` information from the grid all the edges of the cells we have cut out two steps before. We also sort the indices of the edges and delete all duplicates via `np.unique()`. We subtract `1` to get back to python thinking." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([6427302, 6427311, 6427312, 6427313, 6427314, 6427315, 6427316,\n", " 6427317, 6427318, 6427352, 6427353, 6427354, 6427355, 6427356,\n", " 6427357, 6427358, 6427359, 6427360, 6427361, 6427362, 6427363,\n", " 6427365, 6427366, 6427367, 6427368, 6427369, 6427370, 6427371,\n", " 6427372, 6427373, 6427374, 6427375, 6427376, 6427377, 6427385,\n", " 6427387, 6427388, 6427389, 6427390, 6427425, 6427426, 6427482,\n", " 6427483, 6427484, 6427485, 6427486, 6427487, 6427488, 6427489,\n", " 6427490, 6427491, 6427492, 6427493, 6427494, 6427495, 6427496,\n", " 6427497, 6427498, 6427499, 6427500, 6427501, 6427504, 6427507,\n", " 6427508, 6427513, 6427516, 6427527, 6427528, 6427529, 6427531,\n", " 6427532, 6427533, 6427534, 6427535, 6427536, 6427537, 6427538,\n", " 6427539, 6427540, 6427541, 6427542, 6427543, 6427544, 6427545,\n", " 6427546, 6427547, 6427548, 6427549, 6427550, 6427553, 6427602,\n", " 6428151, 6428152, 6428160, 6428161, 6428162, 6428164, 6428165,\n", " 6428166, 6428167, 6428170, 6428202, 6428203, 6428204, 6428205,\n", " 6428206, 6428207, 6428208, 6428213, 6428214, 6428215, 6428216,\n", " 6428217, 6428218, 6428219, 6428410, 6428418, 6428419, 6428420],\n", " dtype=int32)" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "window_edge_indices = (\n", " np.unique(grid.edge_of_cell.isel(cell=window_cell_indices).values) - 1\n", ")\n", "window_edge_indices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Constructing New Grid with Selected Cells, Vertices and Edges" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Wow, that's already great ! We have received a lot of information in the form of indices in individual arrays about our green window. We merge them into one dataset so that everything is compact:" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
       "Dimensions:  (cell: 70, vertex: 50, edge: 119)\n",
       "Coordinates:\n",
       "  * cell     (cell) int64 4282376 4282377 4282378 ... 4283128 4283129 4283130\n",
       "  * vertex   (vertex) int32 2144929 2144932 2144933 ... 2145238 2145239 2145286\n",
       "  * edge     (edge) int32 6427302 6427311 6427312 ... 6428418 6428419 6428420\n",
       "Data variables:\n",
       "    *empty*
" ], "text/plain": [ "\n", "Dimensions: (cell: 70, vertex: 50, edge: 119)\n", "Coordinates:\n", " * cell (cell) int64 4282376 4282377 4282378 ... 4283128 4283129 4283130\n", " * vertex (vertex) int32 2144929 2144932 2144933 ... 2145238 2145239 2145286\n", " * edge (edge) int32 6427302 6427311 6427312 ... 6428418 6428419 6428420\n", "Data variables:\n", " *empty*" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "selected_indices = xr.Dataset(\n", " {\n", " \"cell\": (\"cell\", window_cell_indices),\n", " \"vertex\": (\"vertex\", window_vertex_indices),\n", " \"edge\": (\"edge\", window_edge_indices),\n", " }\n", ")\n", "\n", "selected_indices" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It could be that we need more variables for future calculations. Therefore we create a dictionary with further interesting variables, which we reindex as a precaution." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "vars_to_renumber = {\n", " \"cell\": [\n", " \"adjacent_cell_of_edge\",\n", " \"cells_of_vertex\",\n", " \"neighbor_cell_index\",\n", " \"cell_index\",\n", " ],\n", " \"vertex\": [\"vertex_of_cell\", \"edge_vertices\", \"vertices_of_vertex\"],\n", " \"edge\": [\"edge_of_cell\", \"edges_of_vertex\"],\n", "}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**We now come to the heart of this script: the reindexing.**
\n", "Several things happen here, which is why it is best to define a function. This function `reindex_grid()` needs 3 inputs and returns 1 output. The inputs are the original, complete grid and the parts of the grid that should be reindexed, hence `indices` and `vars_to_renumber`. The `indices` define the cells, vertices and edges. The `vars_to_renumber` are all variables that we are still interested in and can be composed of cells, vertices and edges. Output of our function will be a `new_grid` containing all indices and variables for our green window around the birthplace of Max Planck in such a way that everything starts counting at `0`." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's go through it step by step:
\n", "\n", "**Line 1**: We define a function wiht 3 input variables.\n", "\n", "**Line 2**: We define as `new_grid` the area in the old grid that contains the indices we selected at the beginning of this script for the cells, vertices and edges. For this we use the `.load()` function, which loads the 17GB file into memory and processes it there: this is a little faster.\n", "\n", "**Line 3**: We open a for-loop that accesses the coordinates and the entries of the array `selected_indices`.\n", "\n", "**Line 4**: We open an array, which is only filled with `-2` (exceptional value like `nan` but as an integer) in the original, old grid length and call it `renumbering`.\n", "\n", "**Line 5**: We start counting at `0` at the index positions of the long renumbering array, which belong to the indices of the selected dark green area, until we have reached the length of the short, previously selected array. So what we get is an array with the length of the original grid dimension (20971520 cells, 10485762 vertices, 31457280 edges), which contains a value other than `-2` only at that position within the array which is inside the dark green selected window.\n", "\n", "**Line 6**: we open another for loop over the remaining variables `vars_to_renumber` to be reindexed.\n", "\n", "**Line 7**: For the variables stored in the dictionary of a particular dimension (`cell`, `vertex`, `edge`), we take one item and access it in new_grid (line 2) and subtract `1` to work in python 0-based system; this is done in the square brackets on the right side of the equal sign. We use this to select the valid position in the `renumbering` array but in total we add `1` to output the `new_grid` in the same 1-based thinking as the original `grid`.\n", "\n", "**Line 8**: We output the `new_grid`." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def reindex_grid(grid, indices, vars_to_renumber):\n", " new_grid = grid.load().isel(\n", " cell=indices.cell, vertex=indices.vertex, edge=indices.edge\n", " )\n", " for dim, idx in indices.coords.items():\n", " renumbering = np.full(grid.dims[dim], -2, dtype=\"int\")\n", " renumbering[idx] = np.arange(len(idx))\n", " for name in vars_to_renumber[dim]:\n", " new_grid[name].data = renumbering[new_grid[name].data - 1] + 1\n", " return new_grid" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After long theory we want to use our function and create the actual `new_grid`:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.Dataset>\n",
       "Dimensions:                         (cell: 70, nv: 3, vertex: 50, ne: 6,\n",
       "                                     edge: 119, no: 4, nc: 2,\n",
       "                                     max_stored_decompositions: 4, two_grf: 2,\n",
       "                                     cell_grf: 14, max_chdom: 1, edge_grf: 24,\n",
       "                                     vert_grf: 13)\n",
       "Coordinates:\n",
       "    clon                            (cell) float64 0.1804 0.1805 ... 0.174\n",
       "    clat                            (cell) float64 0.9453 0.946 ... 0.9488\n",
       "    vlon                            (vertex) float64 0.1815 0.1807 ... 0.1717\n",
       "    vlat                            (vertex) float64 0.9456 0.9466 ... 0.9482\n",
       "    elon                            (edge) float64 0.1811 0.1814 ... 0.1724\n",
       "    elat                            (edge) float64 0.9461 0.9471 ... 0.9487\n",
       "  * cell                            (cell) int64 4282376 4282377 ... 4283130\n",
       "  * vertex                          (vertex) int32 2144929 2144932 ... 2145286\n",
       "  * edge                            (edge) int32 6427302 6427311 ... 6428420\n",
       "Dimensions without coordinates: nv, ne, no, nc, max_stored_decompositions,\n",
       "                                two_grf, cell_grf, max_chdom, edge_grf, vert_grf\n",
       "Data variables: (12/91)\n",
       "    clon_vertices                   (cell, nv) float64 0.1794 0.1802 ... 0.173\n",
       "    clat_vertices                   (cell, nv) float64 0.9457 0.9446 ... 0.9492\n",
       "    vlon_vertices                   (vertex, ne) float64 9.969e+36 ... 9.969e+36\n",
       "    vlat_vertices                   (vertex, ne) float64 9.969e+36 ... 9.969e+36\n",
       "    elon_vertices                   (edge, no) float64 0.1807 0.1805 ... 0.1728\n",
       "    elat_vertices                   (edge, no) float64 0.9466 0.946 ... 0.9485\n",
       "    ...                              ...\n",
       "    edge_dual_normal_cartesian_x    (edge) float64 -0.6706 -0.7373 ... -0.7362\n",
       "    edge_dual_normal_cartesian_y    (edge) float64 -0.5071 0.4974 ... 0.4973\n",
       "    edge_dual_normal_cartesian_z    (edge) float64 0.5414 0.4572 ... 0.4589\n",
       "    cell_circumcenter_cartesian_x   (cell) float64 0.576 0.5755 ... 0.5739\n",
       "    cell_circumcenter_cartesian_y   (cell) float64 0.105 0.105 ... 0.1003 0.1009\n",
       "    cell_circumcenter_cartesian_z   (cell) float64 0.8107 0.8111 ... 0.8127\n",
       "Attributes: (12/43)\n",
       "    title:                    ICON grid description\n",
       "    institution:              Max Planck Institute for Meteorology/Deutscher ...\n",
       "    source:                   git@git.mpimet.mpg.de:GridGenerator.git\n",
       "    revision:                 d00fcac1f61fa16c686bfe51d1d8eddd09296cb5\n",
       "    date:                     20180529 at 222250\n",
       "    user_name:                Rene Redler (m300083)\n",
       "    ...                       ...\n",
       "    topography:               modified SRTM30\n",
       "    subcentre:                1\n",
       "    number_of_grid_used:      15\n",
       "    history:                  Thu Aug 16 11:05:44 2018: ncatted -O -a ICON_gr...\n",
       "    ICON_grid_file_uri:       http://icon-downloads.mpimet.mpg.de/grids/publi...\n",
       "    NCO:                      netCDF Operators version 4.7.5 (Homepage = http...
" ], "text/plain": [ "\n", "Dimensions: (cell: 70, nv: 3, vertex: 50, ne: 6,\n", " edge: 119, no: 4, nc: 2,\n", " max_stored_decompositions: 4, two_grf: 2,\n", " cell_grf: 14, max_chdom: 1, edge_grf: 24,\n", " vert_grf: 13)\n", "Coordinates:\n", " clon (cell) float64 0.1804 0.1805 ... 0.174\n", " clat (cell) float64 0.9453 0.946 ... 0.9488\n", " vlon (vertex) float64 0.1815 0.1807 ... 0.1717\n", " vlat (vertex) float64 0.9456 0.9466 ... 0.9482\n", " elon (edge) float64 0.1811 0.1814 ... 0.1724\n", " elat (edge) float64 0.9461 0.9471 ... 0.9487\n", " * cell (cell) int64 4282376 4282377 ... 4283130\n", " * vertex (vertex) int32 2144929 2144932 ... 2145286\n", " * edge (edge) int32 6427302 6427311 ... 6428420\n", "Dimensions without coordinates: nv, ne, no, nc, max_stored_decompositions,\n", " two_grf, cell_grf, max_chdom, edge_grf, vert_grf\n", "Data variables: (12/91)\n", " clon_vertices (cell, nv) float64 0.1794 0.1802 ... 0.173\n", " clat_vertices (cell, nv) float64 0.9457 0.9446 ... 0.9492\n", " vlon_vertices (vertex, ne) float64 9.969e+36 ... 9.969e+36\n", " vlat_vertices (vertex, ne) float64 9.969e+36 ... 9.969e+36\n", " elon_vertices (edge, no) float64 0.1807 0.1805 ... 0.1728\n", " elat_vertices (edge, no) float64 0.9466 0.946 ... 0.9485\n", " ... ...\n", " edge_dual_normal_cartesian_x (edge) float64 -0.6706 -0.7373 ... -0.7362\n", " edge_dual_normal_cartesian_y (edge) float64 -0.5071 0.4974 ... 0.4973\n", " edge_dual_normal_cartesian_z (edge) float64 0.5414 0.4572 ... 0.4589\n", " cell_circumcenter_cartesian_x (cell) float64 0.576 0.5755 ... 0.5739\n", " cell_circumcenter_cartesian_y (cell) float64 0.105 0.105 ... 0.1003 0.1009\n", " cell_circumcenter_cartesian_z (cell) float64 0.8107 0.8111 ... 0.8127\n", "Attributes: (12/43)\n", " title: ICON grid description\n", " institution: Max Planck Institute for Meteorology/Deutscher ...\n", " source: git@git.mpimet.mpg.de:GridGenerator.git\n", " revision: d00fcac1f61fa16c686bfe51d1d8eddd09296cb5\n", " date: 20180529 at 222250\n", " user_name: Rene Redler (m300083)\n", " ... ...\n", " topography: modified SRTM30\n", " subcentre: 1\n", " number_of_grid_used: 15\n", " history: Thu Aug 16 11:05:44 2018: ncatted -O -a ICON_gr...\n", " ICON_grid_file_uri: http://icon-downloads.mpimet.mpg.de/grids/publi...\n", " NCO: netCDF Operators version 4.7.5 (Homepage = http..." ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "new_grid = reindex_grid(grid, selected_indices, vars_to_renumber)\n", "new_grid" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let's see if everything worked as we wanted it and choose reindexed variable like `.vertex_of_cell` and sort it:" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "array([ 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17,\n", " 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34,\n", " 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50])" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "np.unique(new_grid.vertex_of_cell)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Voilà ! That worked and we are done and have now built a new grid-file tailored to the area of our interest, which was provided with a new indexing." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For further processing, the two datasets `selected_indices` and `new_grid` can be saved to have them quickly accessible for further calculations." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "selected_indices.to_netcdf(\n", " f\"selected_indices_region_{bottom_bound}-{top_bound}_{left_bound}-{right_bound}.nc\",\n", " mode=\"w\",\n", ")\n", "new_grid.to_netcdf(\n", " f\"new_grid_region_{bottom_bound}-{top_bound}_{left_bound}-{right_bound}.nc\",\n", " mode=\"w\",\n", ")" ] } ], "metadata": { "jupytext": { "formats": "ipynb,md:myst" }, "kernelspec": { "display_name": "Python 3 (based on the module python3/2022.01)", "language": "python", "name": "python3_2022_01" }, "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.9.9" } }, "nbformat": 4, "nbformat_minor": 4 }