{ "cells": [ { "cell_type": "markdown", "id": "27fede2a-cd34-46b9-9b06-776112f2e680", "metadata": { "tags": [] }, "source": [ "# Processing with HDF5 datasets\n", "\n", "Unlike the `h5py` package, which returns `numpy.ndarray` when accessing the values of datasets, the `h5rdmtoolbox` returns [`xarray.DataArray` objects](https://xarray.pydata.org/). The `xarray.DataArray` object allows to carry attributes with the numpy-like multi-dimensional array. It also supports the concept of dimensions and coordinates, allowing to assign the array axis with meaning ful (meta) data.\n", "\n", "Let's dive into it and explore the practical implications of retrieving `xarray.DatArray`:" ] }, { "cell_type": "code", "execution_count": 1, "id": "9fc2b312-16d5-4aa7-8f10-9107a2e53c08", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "using(\"h5py\")" ] }, "execution_count": 1, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import h5rdmtoolbox as h5tbx\n", "import numpy as np\n", "\n", "h5tbx.use(None)" ] }, { "cell_type": "markdown", "id": "8bc428f1-f2a6-4f4d-bbbf-9880358cff06", "metadata": {}, "source": [ "Let's create an example file. Note, that we pass `make_scale` and `attach_scale` as arguments to setup the coordinates and their association to the HDF5 dataset \"data\". The useful implications will be visible when we access the dataset values in the next steps." ] }, { "cell_type": "code", "execution_count": 2, "id": "bf629a23-0162-4d0b-a910-691b4ec97941", "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "

\n", "\n", "

\n", " \n", " \n", " /(3)\n", " \n", "
\n", "\n", "
\n", "\n", "
\n", "\n", "
\n", "

\n", "

" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "with h5tbx.File() as h5:\n", " dsx = h5.create_dataset('x', data=np.linspace(0, 10, 5),\n", " attrs=dict(units='mm', long_name='x'),\n", " make_scale=True)\n", " dsy = h5.create_dataset('y', data=np.linspace(0, 5, 11),\n", " attrs=dict(units='mm', long_name='y'),\n", " make_scale=True)\n", " h5.create_dataset('vel', data=np.random.random((11, 5)),\n", " attrs=dict(units='m/s', long_name='velocity'),\n", " attach_scales=(dsy, dsx))\n", " h5.dump()" ] }, { "cell_type": "markdown", "id": "e2b2da86-50d5-48c2-a15a-7fce8e06fc66", "metadata": {}, "source": [ "## Array Slicing\n", "\n", "Slicing an HDF5 dataset returns a `xarray.DataArray`" ] }, { "cell_type": "code", "execution_count": 3, "id": "a9bd4492-84ad-4950-81ea-eae066c591e8", "metadata": {}, "outputs": [ { "data": { "text/html": [ "

\n", "

<xarray.DataArray 'vel' (y: 11, x: 5)>\n",
       "0.5832 0.7347 0.3899 0.7013 0.07984 ... 0.5462 0.07077 0.3864 0.2998 0.2473\n",
       "Coordinates:\n",
       "  * y        (y) float64 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0\n",
       "  * x        (x) float64 0.0 2.5 5.0 7.5 10.0\n",
       "Attributes:\n",
       "    long_name:  velocity\n",
       "    units:      m/s

" ], "text/plain": [ "\n", "0.5832 0.7347 0.3899 0.7013 0.07984 ... 0.5462 0.07077 0.3864 0.2998 0.2473\n", "Coordinates:\n", " * y (y) float64 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0\n", " * x (x) float64 0.0 2.5 5.0 7.5 10.0\n", "Attributes:\n", " long_name: velocity\n", " units: m/s" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "with h5tbx.File(h5.hdf_filename) as h5:\n", " data = h5.vel[:]\n", "data" ] }, { "cell_type": "markdown", "id": "87e1093b-3e45-4d01-8ab6-42a58a2cdf6e", "metadata": {}, "source": [ "## Advantages of retrieving `xarray.DataArray`\n", "\n", "Few of the advantages:\n", "- attributes (aux. info) with the array `.attrs` (copied from the HDF dataset to the `xarray.DataArray`)\n", "- dimensions and coordinates (1D arrays) to address the axis by label rather than by idex\n", "- apply operations (computations, visualizations) based on the meta data" ] }, { "cell_type": "code", "execution_count": 4, "id": "745ae071-ef36-46b8-89ec-4df1bc02da9e", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'long_name': 'velocity', 'units': 'm/s'}" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.attrs" ] }, { "cell_type": "markdown", "id": "4dbc2472-9494-412d-9d16-396ca8eb2dba", "metadata": {}, "source": [ "Select subarray by specifying coordinate values for a given axis (coordinate):" ] }, { "cell_type": "code", "execution_count": 5, "id": "a05311e4-357e-4832-8cf5-752a6de17241", "metadata": {}, "outputs": [ { "data": { "text/html": [ "

\n", "

<xarray.DataArray 'vel' (x: 5)>\n",
       "0.3443 0.3647 0.09872 0.8726 0.3015\n",
       "Coordinates:\n",
       "    y        float64 2.0\n",
       "  * x        (x) float64 0.0 2.5 5.0 7.5 10.0\n",
       "Attributes:\n",
       "    long_name:  velocity\n",
       "    units:      m/s

" ], "text/plain": [ "\n", "0.3443 0.3647 0.09872 0.8726 0.3015\n", "Coordinates:\n", " y float64 2.0\n", " * x (x) float64 0.0 2.5 5.0 7.5 10.0\n", "Attributes:\n", " long_name: velocity\n", " units: m/s" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.sel(y=2.0)" ] }, { "cell_type": "markdown", "id": "fb30dc33-8517-445f-88b9-ac4b9118e339", "metadata": {}, "source": [ "Plot data by using information from attributes and coorinates:" ] }, { "cell_type": "code", "execution_count": 6, "id": "a4799387-5349-4083-9e79-aa74bb639fbe", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiUAAAG2CAYAAACgd/abAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABFG0lEQVR4nO3de1hUdeLH8c+AclkRFBUQBbEszbyVKOEla3XXzCy2Z1szU7zUbuWdrVVcBd1StNIl05Uuv6x219Vulrql66IiPmteo7XHn5qpaRqolaAYoMz8/vDHxDADzJw5c77f7zmf1/PwPDLOnPnCnDnnzfecmbE5HA4HiIiIiAQLEj0AIiIiIoBRQkRERJJglBAREZEUGCVEREQkBUYJERERSYFRQkRERFJglBAREZEUGCVEREQkBUYJERERSYFRQkRERFIQGiXz5s2DzWZz+erSpYvIIRERESlpx44dGDFiBOLj42Gz2fDhhx82epvt27fj9ttvR2hoKDp16oQ333wz4ONsiPCZkltvvRXffvut82vnzp2ih0RERKSc8vJy9OzZEytWrPDq+idOnMDw4cNx9913o6ioCNOnT8djjz2GzZs3B3ik9Wsi7J5rBtCkCeLi4kQPg4iISGnDhg3DsGHDvL5+Xl4eOnbsiCVLlgAAbrnlFuzcuRN//vOfMXTo0EANs0HCo+TLL79EfHw8wsLCkJqaipycHCQmJnq8bmVlJSorK53f2+12fP/992jVqhVsNptRQyYiIgU5HA5cunQJ8fHxCAoK3IGCiooKVFVV+b0ch8Phtm8LDQ1FaGio38sGgF27dmHIkCEulw0dOhTTp0/XZflaCI2SlJQUvPnmm+jcuTO+/fZbzJ8/HwMHDsQXX3yB5s2bu10/JycH8+fPFzBSIiIyi9OnT6N9+/YBWXZFRQU6dohA8blqv5cVERGBy5cvu1yWnZ2NefPm+b1sACguLkZsbKzLZbGxsSgrK8OPP/6I8PBwXe7HF0KjpPY0U48ePZCSkoIOHTrgnXfewcSJE92un5mZiYyMDOf3paWlSExMRML8OQgKCzNkzL6I2+kQPYRGhW/YJ3oIfvlxRLLoIRCRIq5drcD+TQs8/tGrl6qqKhSfq8bX+5MQ2Vz7bEzZJTs69D6J06dPIzIy0nm5XrMkshJ++Ka2Fi1a4Oabb8axY8c8/n9901ZBYWFSRsk511kxxBfIFylXHxzg8n34ut2CRtKwH3+V4vFyqVZgwc4O8u0QpozrI5ERjDjcH9k8CJHNg/1fTmSkS5ToKS4uDiUlJS6XlZSUIDIyUsgsCSDZNv3y5cv46quvMGbMGNFDCYi6Ow0Zdwp1d/5GR0p98UHX+RoegVyWjOsvEXkvNTUVH3/8sctlW7ZsQWpqqqARCY6Sp59+GiNGjECHDh1w9uxZZGdnIzg4GKNGjRI5LMNYOVIYH/XTMzwCyZ9xyriuE6nu8uXLLkcaTpw4gaKiIkRHRyMxMRGZmZk4c+YM3n77bQDAE088geXLl+MPf/gDJkyYgK1bt+Kdd97BP//5T1E/gtgo+eabbzBq1Ch89913aNOmDQYMGIBPP/0Ubdq0ETksYcwYKYyP+qkSH4HAWRoi/e3btw9333238/uaczDT09Px5ptv4ttvv8WpU6ec/9+xY0f885//xIwZM/DSSy+hffv2eP3114W9HBgAbA6HQ9lnd1lZGaKiotBh8XNSnlOiN26I1WPl8JAZn0vWdO1qBXZvmIvS0tKAnadRs1/64egNfp1TUnapGi1vPh7QscpIqnNKqGEqzKRYEcNDPZypIZKTKaIkNfkwQiJCnN8XftpV4GiMw0gxFuODamhdF/gcJWqYKaKkroF3HHK7zAqhwkjxH8ODAqm+9YvPVaLrTBklnlgxVBgpnjE8SDYNrZN83pKVWCZKPLFaqFgpUhgeZBacXSErsXSUeGKlUDFDpDA+yKo4u0JmxCjxglVCRdZIYXgQ+YazK6QqRolGnkIFMFesGBkpDA+iwOPsCsmOUaIzM8+q+BspDA8ieXF2hWTAKDGAWUOlvkhhfBCZB2dXyEiMEkHMGCqMESJrYbCQ3hglEjFjqBCRNfFwEGnBKJEcQ4WIzISzK9QQRomCGCpEZEacXSFGiUkwVIjIrDi7Yh2MEhNjqBCR2XF2xVwYJRbDUCEiK/AUK/YKG7BBwGDIa4wSYqgQEZEUGCXkEUOFiIiMxighrzFUiNzV9zlYjeFzh8gdo4T80tAGmRtdUp3W4NCybD5fiBglFEC+bNC5QSajBTI4tOBMJBGjhCTh7Q6CG2nyhmzBoRVDhayGUUJK4ewLmSU4tGKokJkxSsi0OPuiJqtHhxY8P4XMglFClsd4MQ6DwxicTSFVMUqIvMRDRw1jcMiNoUIqYJQQBYDZZl8YHObEwz4kG0YJkUAyzL4wOKgGZ1NINEYJkSK0zL4wOMhfDBUyUpDoARCRvhgiRKQqzpQQKcTXv1C9uT4jhhrDmREyCqOESAIiN/oMF6rB+CDRGCVEAWSWjby3PwfjRR1mWTfJXBglRD7ixrx+jBc5cZ0lVTBKiP4fN9zG4SGjwOF6TCpjlJDpcSOtJs66NI7rNpkNo4SUxI0x1bDKrAvXebICRglJhRteCgTVZl34PCCrYpSQIbiRJRUYHS98XhC5YpSQbriBJavQcsiIzw+ixjFKyC/c0BJ5xucGke8YJeQzbmyJiCgQGCXkNcYIEREFEqOEGsQQISIiozBKyA1DhIiIRGCUkBNjhIiIRGKUWBxDhIiIZMEosSjGCBERyYZRYiEMESIikhmjxOQYIkREpApGiUkZGSPxBQ6Pl58dZDNsDEREpD5GiYkYPStSX4w09P8MFSIiqg+jRHGyhYivt2ekEBFRDUaJomQ4PBOIZTNSiIisi1GiENVmRfS6T4YKEZE1MEokZ4UQaQxnU4iIrIFRIimzHJ4JBEYKycyb5xPXWSLPGCUS4ayINjzkQ6JofQ4xrIk8Y5RIgLMi+uNGn/TGE76JAo9RIghDxFjc6JMvRD9nuL6SVTFKDGTmwzPh63Y7//3jr1IMu1+teMiHaogOEG8wUsgqGCUGMPOsSO0YaegyFUOFG37zUSFAvMF1lcwqSPQAaixatAg2mw3Tp08XPRRdFH7a1fllhPgCh2Eb3PB1u51fWm7jy+1EqvmdGvm7Jf3UffzM/Bha6Welhq1YsQJJSUkICwtDSkoK9uzZ0+D1c3Nz0blzZ4SHhyMhIQEzZsxARUWFQaN1J8VMyd69e/HKK6+gR48eoofiF6scngnE8lScSQH4F6osuBN2x9kU61m7di0yMjKQl5eHlJQU5ObmYujQoThy5AhiYmLcrr969WrMmjULb7zxBvr164ejR49i3LhxsNlsWLp0qYCfQIIouXz5MkaPHo3XXnsNzz33nOjhaGK1wzNG3Y+KocINf+AxQLThump+S5cuxeOPP47x48cDAPLy8vDPf/4Tb7zxBmbNmuV2/f/85z/o378/HnnkEQBAUlISRo0ahd27xc1mC4+SSZMmYfjw4RgyZEijUVJZWYnKykrn92VlZYEeXr04KxJ4ZphN4YbfPwyQwOG6qoa6+7nQ0FCEhoa6Xa+qqgr79+9HZmam87KgoCAMGTIEu3bt8rjsfv364W9/+xv27NmDvn374vjx4/j4448xZswYfX8IHwiNkjVr1uDAgQPYu3evV9fPycnB/PnzAzyq+pk5RGrIEiSeqPYKH4CHfHzBABGr5vfP9VMfT57uh5CIEM23r7pcBeA4EhISXC7Pzs7GvHnz3K5/4cIFVFdXIzY21uXy2NhYHD582ON9PPLII7hw4QIGDBgAh8OBa9eu4YknnsDs2bM1j9tfwqLk9OnTmDZtGrZs2YKwsDCvbpOZmYmMjAzn92VlZW4PWKCY+RBNbT/+KkXqMFElRhoSX+Dghv//MUTkwPVRXqdPn0ZkZKTze0+zJFpt374dCxcuxF/+8hekpKTg2LFjmDZtGp599lnMnTtXt/vxhbAo2b9/P86dO4fbb7/deVl1dTV27NiB5cuXo7KyEsHBwS63qW/aKpCsEiO1yRgmZoiR2qz+V6ks67rVWXX9U0lkZKRLlNSndevWCA4ORklJicvlJSUliIuL83ibuXPnYsyYMXjssccAAN27d0d5eTl++9vf4o9//COCgox/ga6wKBk8eDAOHjzoctn48ePRpUsXzJw50y1IjGbFGKlNljAxW4zUZaU4kXE9tyIrrGtWFBISgt69eyM/Px9paWkAALvdjvz8fEyePNnjba5cueIWHjX7XodDzPNVWJQ0b94c3bp1c7msWbNmaNWqldvlRrJ6jNQmOkzMHiS1mfWQjuzruFWYcd0idxkZGUhPT0dycjL69u2L3NxclJeXO1+NM3bsWLRr1w45OTkAgBEjRmDp0qW47bbbnIdv5s6dixEjRgibGBD+6htZWOEkVi1EhImVYqQ2s8yaqLJum53q6xH5buTIkTh//jyysrJQXFyMXr16YdOmTc6TX0+dOuUyMzJnzhzYbDbMmTMHZ86cQZs2bTBixAgsWLBA1I8Am0PUHI0OysrKEBUVhYfzH9V8ljNjxDtGhIlVY8QTFXcoqq7bZqPiumMUe0UFvp45B6WlpV6dp6GFHvsl4Pqrb9YM/ltAxyojy86UMEZ8E8gZE8aIO1VmTVRfr81C9vWEyFuWixLGiHZ6hwljpHEynmtipnVaZbKtF0R6sEyUMEb0oVeYMEi8J8OsiVnXZxUxRsjMLBElfEWNvvwJE8aIdiLixArrswoYImQVpo4Sxkjg1MSFt3HCGNFPoA/pWG1dlhVDhKzIlFHCGDFOY7MmjJHA0HvWxOrrsSwYImR1pooSxogYnsKEMWIMf2ZNuA7LgzFCdJ0pomTXvi4I8vJD/fzFDblntcOEQWIsX2dNuA7LgSFC5M4UUWIUbswbxhgRq6FZE667cmCIEDWMUeIFbtBJFbVnTbjeyoEhQuQ9RkkDuFEnVXHdFY8xQuQ7RokH3KATkRYMESL/MEpqYYwQkRaMETXE7XTga9GDoAYxSsAYISLfMUTUUHv7fk3gOMg7lo4SxggR+YIhog5u39VkySjhykpEvmCMqIHbdvVZKkq4whKRtxgiauB23VwsEyVccYmoMQwRdXCbbk6mjxKuuETUGMaIGrg9Nz/TRglXXiJqCENEDdyWW4vpooQrMBHVhyGiDm7Lrck0UcIVmIjqwxhRA7fjZIooidvpAJqKHgURyYQhogaGCNVmiighIqrBGFEDY4Q8YZQQkfIYImpgiFBjGCVEpCSGiDoYI+QtRgkRKYUxogaGCGnBKCEi6TFE1MAQIX8xSohISgwRdTBGSC+MEiKSCmNEDQwRCgRGCREJxxBRT81jxjghPTFKiEgIhog5eHocGSqkFaOEiAzDELGGuo8zI4W8xSghXfCvJaoPQ4S4fSBvMUrILw3tcBrbGXGjZG6MEWoIQ4U8YZSQJnrscBpaBjdOamKIkD942IcYJeQTo3Y6nGVRB0OEAoWzKdbDKCGvybTzYbSIJdO6QNbCUDE3Rgk1SsUdEA8N6U/F9YCsgYd9zINRQvUy606IsyzeM+s6QObG2RR1MUrIjdV3RIwWrgNkPmcH2WCvsAEbRI+EGsIoISfuiLxj1mjh409EojFKiDsjnal0PgsfeyKSCaPEwrhDMp4Msyx83IlIVowSC+JOSV6BihY+5kSkAkaJhXDHpD5fooWPNxGphlFiEdxBWQMfZyJSGaPE5LiTIiIiVTBKTIoxQkREqmGUmAxjhEi8gXccCvh9FH7aNeD3QWQ0RolJMEaIxDIiRIy4P8YOicQoURxjhEgMoyPEKIH8uRg81BhGiaIYI0TGMmuEGImzO9QYRoliGCNExmGIqMHbx6nqchW+DvBYyD+MEkUwRogCjxFCJBajRHKqxIinjTmnVEkFDBEieTBKJKZCkDS0QW9sY89oIREYIUTyYpRISPUY0WMZDBbSCyOESB2MEolYJUb8vR8GCzWGIUKkJkaJBBgjvuFhIapLpvWTyMxuv/12n65vs9mwfv16tGvXzqvrM0oEUiFGAPU2+JxlsQbV1ksiMygqKsLvf/97RERENHpdh8OBRYsWobKy0uvlM0oEYIyIw2BRlxnXRyIVPfPMM4iJifHqukuWLPFp2YwSAzFG5MbDQvKx6rpIJKsTJ06gTZs2Xl//0KFDiI+P9/r6QVoGpZeVK1eiR48eiIyMRGRkJFJTU/HJJ5+IHFJAnB1kUyJIBt5xiDuBBtT8fjx9kT74eyXyz4oVK5CUlISwsDCkpKRgz549DV7/4sWLmDRpEtq2bYvQ0FDcfPPN+Pjjj+u9focOHWCzeb8/S0hIQHBwsNfXFzpT0r59eyxatAg33XQTHA4H3nrrLTzwwAP47LPPcOutt4ocmi5UCBGAf43qgYeFtOG6R6SftWvXIiMjA3l5eUhJSUFubi6GDh2KI0eOeDzcUlVVhV/84heIiYnBe++9h3bt2uHrr79GixYtvLq/TZs2ISIiAgMGDABwPYhee+01dO3aFStWrEDLli19/hlsDofD4fOtAig6OhovvPACJk6c2Oh1y8rKEBUVhZQRz6JJ0zADRucdxgj5wmrRwvWORKm6XIU1g/+G0tJSREZGBuQ+avZLD+c/ipCIEM3L0TLWlJQU9OnTB8uXLwcA2O12JCQkYMqUKZg1a5bb9fPy8vDCCy/g8OHDaNq0qc9j7N69OxYvXox7770XBw8eRJ8+fZCRkYFt27ahS5cuWLVqlc/LlOackurqarz77rsoLy9Hamqqx+tUVla6nMVbVlZm1PC8pkKQcKcgF7PPsnB9I5FWJRY6/112qRprBI5Fi7r7udDQUISGhrpdr6qqCvv370dmZqbzsqCgIAwZMgS7du3yuOz169cjNTUVkyZNwkcffYQ2bdrgkUcewcyZM7065HLixAl07Xp9G/X+++/jvvvuw8KFC3HgwAHce++9vvyYTsKj5ODBg0hNTUVFRQUiIiKwbt065w9ZV05ODubPn2/wCL2nQpCQWswQLIWfdmWYkOXs2tcFQWHaZ/DtFRUArp+TUVt2djbmzZvndv0LFy6guroasbGxLpfHxsbi8OHDHu/j+PHj2Lp1K0aPHo2PP/4Yx44dw1NPPYWrV68iOzu70TGGhITgypUrAIB///vfGDt2LIDrRzy0ThoIj5LOnTujqKgIpaWleO+995Ceno6CggKPYZKZmYmMjAzn92VlZW4PGJFV1LejlzFW9B4TI4e8Nf7UQOe/qy5XATgubjAanD592uXwjadZEq3sdjtiYmLw6quvIjg4GL1798aZM2fwwgsveBUlAwYMQEZGBvr37489e/Zg7dq1AICjR4+iffv2msYkPEpCQkLQqVMnAEDv3r2xd+9evPTSS3jllVfcrlvftBUR/aT2DlvGQNEDI4esoubVqY1p3bo1goODUVJS4nJ5SUkJ4uLiPN6mbdu2aNq0qcuhmltuuQXFxcWoqqpCSEjD58QsX74cTz31FN577z2sXLnS+a6tn3zyCe65555Gx+yJ8Cipy263+/Tub0RUv4F3HDJtmOiJkUOqCwkJQe/evZGfn4+0tDQA1/en+fn5mDx5ssfb9O/fH6tXr4bdbkdQ0PV3CDl69Cjatm3bYJBs3boVgwYNQmJiIjZu3Oj2/3/+8581/xxCoyQzMxPDhg1DYmIiLl26hNWrV2P79u3YvHmzyGERmQrDxHgN/b4ZLBQoGRkZSE9PR3JyMvr27Yvc3FyUl5dj/PjxAICxY8eiXbt2yMnJAQA8+eSTWL58OaZNm4YpU6bgyy+/xMKFCzF16tQG7+exxx7DxYsXcc899yAtLQ3Dhg1D8+bNdfkZhEbJuXPnMHbsWHz77beIiopCjx49sHnzZvziF78QOSwi06nZETJOxKv7GDBSSC8jR47E+fPnkZWVheLiYvTq1QubNm1ynvx66tQp54wIcP0k2s2bN2PGjBno0aMH2rVrh2nTpmHmzJkN3s/x48fx3//+F+vXr8eLL76I9PR0DBgwAPfffz8eeOABJCYmav4ZpHufEl/I9j4lqrz6hhtBa2OYyIvPzcAy8n1KOix+zu9X33w9c05Ax6qXs2fPYv369Vi/fj22bduGzp074/7778f999+P5ORkn5Yl9G3mich43PHJq/DTri5fRCqIj4/HE088gY8//hgXLlzAnDlzcPLkSdxzzz1YuHChT8uS7kRXCizukAjg4RxV8FAPqaZZs2b49a9/jV//+teorq7G999/79PtGSVEFsaTYNXCSCHZ7N27F9u2bcO5c+dgt9udl9tsNixZssSnTxQGGCW6UeV8EqK6GCbqYqSQSAsXLsScOXPQuXNnxMbGunx6sC+fJFwbo4SIGCYmwUghI7300kt44403MG7cON2WySghIgA8z8SMGCkUSEFBQejfv7++y9R1aUSkPO64zIuv7iE9zZgxAytWrNB1mV7NlCxbtsznBY8fP163d3gjImPxcI41cCaF/PH0009j+PDhuPHGG9G1a1c0bdrU5f8/+OADn5fpVZRMnz4d7du3d/nQnoacPn0a9913H6OESGEME+thpJAvpk6dim3btuHuu+9Gq1atNJ/cWpvX55Ts27cPMTExXl2XMUJkDjzPxNpqP+4MFKrrrbfewvvvv4/hw4frtkyvoiQ7OxsRERFeL3T27NmIjo7WPCgikgtnTYizKFRXdHQ0brzxRl2X6dWJrtnZ2fjZz37m9UIzMzPRokULrWMiIglxJ0S18aRZmjdvHrKzs3HlyhXdlsmXBBOR13g4h+rDmRTrWbZsGb766ivExsYiKSnJ7UTXAwcO+LxMn6Pku+++Q1ZWlse3lQXg8/vcE5F6eDiHGsNIMb+0tDTdl+lzlIwZMwbHjh3DxIkT3d5W1qpUeYt5bhRITwwT8gUjxXyys7N1X6bPUVJYWIidO3eiZ8+eug+GiNTCMCGtGCnkic9R0qVLF/z444+BGAsRKYjnmTQsvsDh/Lcqs6oiMFLUEB0djaNHj6J169ZeXT8xMRGFhYXo0KGDV9f3OUr+8pe/YNasWcjKykK3bt3cTmyJjIz0dZFEZAKcNXFVO0bqu4yRUj9GipwuXryITz75BFFRUV5d/7vvvkN1dbXXy/c5Slq0aIGysjL8/Oc/d7nc4XDAZrP5dOdEZC4Mk+s8BUlj12OgNIyRIo/09PSALdvnKBk9ejSaNm2K1atX80RXIgmtSiwEAIw/NVDI/Vs5TLyNEW9uy0hpGN9tVoy6r7jVm89R8sUXX+Czzz5D586dAzEeIvJBTYA09H8i4sRq55n4EyPeLJOB0jDOopiHz1GSnJyM06dPM0qIDNZQgDR2O86aBE4ggqSx+2CkNIyRoi6fo2TKlCmYNm0annnmGXTv3t3tRNcePXroNjgiq9IaIA0tj2GiLyNixJv7ZqA0rmb9s1dUCB4JNcbnKBk5ciQAYMKECc7LbDYbT3Ql0kjvAGnofkSGCWCOwzkiY8QTzqKQmfgcJSdOnAjEOJTFDQD5wqgAaez+OWuijWxB4glnUUhlPkeJt2+AQkTiI6Q+omdNVAsTFWLEE86iUCANGjQIEydOxEMPPYTw8HBdlqnpU4LPnj2LnTt3evxAvqlTp+oyMCLVyBog9WGYeEfVIPGEkUJ6uu222/D0009jypQp+M1vfoOJEyfijjvu8GuZPkfJm2++id/97ncICQlBq1atXN6nxGazMUokxbPP9aVagNSHLxuun5lipD481EP+yM3NxYsvvoj169fjrbfewp133olOnTphwoQJGDNmDGJjY31eps9RMnfuXGRlZSEzMxNBQUE+3yGRaswSIA3hrMlPrBAjnnAWhbRo0qQJHnzwQTz44IM4d+4cXn31VcydOxezZ8/Gvffei6lTp7q9A3yDy/N1AFeuXMHDDz/MICHTskKEeMIwsW6QeMJZFPLFnj17sGrVKqxZswYxMTEYN24czpw5g/vuuw9PPfUUXnzxRa+W43NZTJw4Ee+++67PAyaxZNjgq0LUjpnE4873J2cH2ZxfRJ6cO3cOS5YsQbdu3TBw4ECcP38e//jHP3Dy5EnMnz8fr7/+Ov71r38hLy/P62X6PFOSk5OD++67D5s2bfL45mlLly71dZFKq/lrQoUnbk2Y8PySxvkaJqrOrogMMFlDufZz2SozJypsv0g+7du3x4033ogJEyZg3LhxaNOmjdt1evTogT59+ni9TE1RsnnzZufbzNc90dWqGCfWpmXnLipkGCLeq3k+myVOVNg+kTry8/MxcGDD25PIyEhs27bN62X6HCVLlizBG2+8gXHjxvl6U0tgnJC3jJyNYYj4R7XZExW2P6S+7OxsfPDBB2jRooXL5WVlZUhLS8PWrVt9XqbPURIaGor+/fv7fEdWwzghval0rosZQqQ+MgWKCtsXMq+CggJUVVW5XV5RUYHCQm1/RPkcJdOmTcPLL7+MZcuWabpDq2GckBWYOUIaYuThHRW2IWQN//3vfwEADocDhw4dQnFxsfP/qqursWnTJrRr107Tsn2Okj179mDr1q3YuHEjbr31VrcTXT/44ANNAzE7xgmZjVVDxBM9Z09U2EaQtfXq1Qs2mw02m83je5CEh4fj5Zdf1rRsn6OkRYsWePDBBzXdGTFOSG0MkcZ5O3uiwjaAyJMTJ07A4XDghhtuwJ49e1xedRMSEoKYmBgEBwdrWrbPUbJq1SpNd0SuGCekAkaIdio8t4m0qPlg3rqffacHTR/IR/phnJBsGCJEVJ/169dj2LBhaNq0KdavX9/gde+//36fl+9VlNx+++3Iz89Hy5YtvVrogAEDsHbtWs0nulgR44REYogQkTfS0tJQXFyMmJgYpKWl1Xs9m82G6upqn5fvVZQUFRXh888/R3R0tFcLLSoqQmVlpc+DIcYJGYMRQkRa1D5kI/TwzeDBg+FweHdWuZXf2VUvjBPSG0OEiGTnVZScOHHC5wW3b9/e59uQOxXjhBpnVMDxMSGiQJk6dSo6deqEqVOnuly+fPlyHDt2DLm5uT4v06soqTnTlsRRKU6ocYwFIlLd+++/7/Fk1379+mHRokWaoiRIh3GRgeILHMLf2pqIiOi7775DVFSU2+WRkZG4cOGCpmUyShTFOCEiIpE6deqETZs2uV3+ySef4IYbbtC0TL5PieJ4WIeIiETIyMjA5MmTcf78eefbzefn52PJkiWaDt0AjBLTYJwQEZGRJkyYgMrKSixYsADPPvssACApKQkrV67E2LFjNS3T58M36enp2LFjh6Y7o8DjYR0iIjLKk08+iW+++QYlJSUoKyvD8ePHNQcJoCFKSktLMWTIENx0001YuHAhzpw5o/nOKXAYJ0RyGHjHIb6HD5na+fPnceTIERQVFWk+wbWGz1Hy4Ycf4syZM3jyySexdu1aJCUlYdiwYXjvvfdw9epVvwZD+quJEwYKkbHqxkjN9wwUMovy8nJMmDABbdu2xZ133ok777wTbdu2xcSJE3HlyhVNy9T06ps2bdogIyMDn3/+OXbv3o1OnTphzJgxiI+Px4wZM/Dll19qGgwFFuOEKPC8CQ8GCplBRkYGCgoKsGHDBly8eBEXL17ERx99hIKCAvz+97/XtEy/TnT99ttvsWXLFmzZsgXBwcG49957cfDgQXTt2hXPP/88ZsyY4c/iKUB4Uqx51Y5OPr7G8Scuat+Wb6oXWHE7Hfha9CBM5P3338d7772Hu+66y3nZvffei/DwcPzmN7/BypUrfV6mz1Fy9epVrF+/HqtWrcK//vUv9OjRA9OnT8cjjzyCyMhIAMC6deswYcIERonkGCfq0DLDxUAJPL1nOmqWxzjRX3yBA9dED8Jkrly5gtjYWLfLY2JiNB++8TlK2rZtC7vdjlGjRmHPnj3o1auX23XuvvtutGjRQtOAyHiME7GMOKTGx1hfgT7swtkTffGwdWCkpqYiOzsbb7/9NsLCwgAAP/74I+bPn4/U1FRNy/Q5Sv785z/joYcecg7AkxYtWmj6ED8Sizsufcm4IeTsiX9EnAPCQNFOxuegmbz00ksYOnQo2rdvj549ewIAPv/8c4SFhWHz5s2alulzlIwZM0bTHZE6GCeNM8PGjoHiPVlOSGWgeM8Mz1HZdevWDV9++SX+/ve/4/DhwwCAUaNGYfTo0QgPD9e0TL6jK9XLqnFixY0ZA8UzWWLEEwZK/az4HBblZz/7GR5//HHdlscooUaZJU64ofKOWR5vrWQOkfowUH7C53lgrV+/3uvr3n///T4vn1FCXpN5Z8UNkf6sNnuiYox4YuVA4XYg8NLS0ry6ns1mQ3V1tc/LZ5SQz4yME25k5GDmQDFLjHhipUDhtsIYdrs9oMtnlJBm/sYJNyJqMkugGBUjqxILAQDjTw005P7qY+ZA4bZEvIqKigZflestoVGSk5ODDz74AIcPH0Z4eDj69euHxYsXo3PnziKHRT6qGyfcQBgnfN3uRq/z469SAnb/Mh/Sq48RMVITIvVdxkDRD7c3rlasWIEXXngBxcXF6NmzJ15++WX07du30dutWbMGo0aNwgMPPIAPP/zQq/uqrq7GwoULkZeXh5KSEhw9ehQ33HAD5s6di6SkJEycONHn8QuNkoKCAkyaNAl9+vTBtWvXMHv2bPzyl7/EoUOH0KxZM5FDIw24cdCfN9Hh7TKMiBNA3kARFSONXY+Bog23N+7Wrl2LjIwM5OXlISUlBbm5uRg6dCiOHDmCmJiYem938uRJPP300xg40Ld1ccGCBXjrrbfw/PPPu7wCp1u3bsjNzdUUJTaHwyHNI3v+/HnExMSgoKAAd955Z6PXLysrQ1RUFFJGPIsmTf2fNiIymh7RoUUgA6U20YFi9CEaf4kOlNpkDhStQXLtagV2b5iL0tJS58ei6K1mv9Rh8XMI8uNwhr2iAl/PnOPTWFNSUtCnTx8sX778+jLsdiQkJGDKlCmYNWuWx9tUV1fjzjvvxIQJE1BYWIiLFy96PVPSqVMnvPLKKxg8eDCaN2+Ozz//HDfccAMOHz6M1NRU/PDDD14tpzapzikpLS0FAERHR3v8/8rKSlRWVjq/LysrM2RcRFqJio7G1B6XGQ/vqBYjnpYnOlBknUGx2gxJ3f1caGgoQkND3a5XVVWF/fv3IzMz03lZUFAQhgwZgl27dtW7/D/96U+IiYnBxIkTUVjo2/p85swZdOrUye1yu92Oq1ev+rSsGtJEid1ux/Tp09G/f39069bN43VycnIwf/58g0dG5BtZQ6Q+4et2B3zmxMg4UTVIPC1fdJiQdnE7HWjSVHtAXbt6/RONExISXC7Pzs7GvHnz3K5/4cIFVFdXu31AXmxsrPPdVuvauXMn/ud//gdFRUWaxti1a1cUFhaiQ4cOLpe/9957uO222zQtU5oomTRpEr744gvs3Lmz3utkZmYiIyPD+X1ZWZnbA0YkWu0dvAqBYsShHNGHccg/A+84JM1sydlBNkvNlpw+fdrl8I2nWRItLl26hDFjxuC1115D69atNS0jKysL6enpOHPmDOx2Oz744AMcOXIEb7/9NjZu3KhpmVJEyeTJk7Fx40bs2LED7du3r/d69U1bEclK9kAxW5CYZZak9v3IMlvCMBEjMjLSq3NKWrdujeDgYJSUlLhcXlJSgri4OLfrf/XVVzh58iRGjBjhvKzmPUiaNGmCI0eO4MYbb2zwPh944AFs2LABf/rTn9CsWTNkZWXh9ttvx4YNG/CLX/zCmx/PjdAocTgcmDJlCtatW4ft27ejY8eOIodDFFCyBQqDRBujgkRGDBN5hYSEoHfv3sjPz3e+66rdbkd+fj4mT57sdv0uXbrg4MGDLpfNmTMHly5dwksvveT1UYiBAwdiy5Ytfo+/htAomTRpElavXo2PPvoIzZs3R3FxMQAgKipK8ycMEqlAdKAwSNQh02wJyS0jIwPp6elITk5G3759kZubi/LycowfPx4AMHbsWLRr1w45OTkICwtzO3+zRYsWAFDveZ11PfbYY3j00Udx11136fYzCI2SlStXAoDbD7Rq1SqMGzfO+AERCVA3EAIZKWZ8KbCRQSJqlkSmMOFsibxGjhyJ8+fPIysrC8XFxejVqxc2bdrkPPn11KlTCAoK0u3+zp8/j3vuuQdt2rTBww8/jNGjR6NXr15+LVOq9ynxFd+nhMxOz0BhkPhH9GEbWaKkhixhAnj/MmEj36fE3/2SEWPVww8//IB3330Xq1evRmFhIbp06YLRo0fjkUceQVJSks/L0y+ZiEh3P/4qxfnl73KMYNYgkYHoKKpLpt8/X90lTsuWLfHb3/4W27dvx9dff41x48bhr3/9q8f3L/EGo4R0ww1DYGkNFAaJ/2QJAlnGISNuf8S6evUq9u3bh927d+PkyZNu75fiLSleEkxqq70xqLth4PHewPD2RFkGCQWSTOeXkBjbtm3D6tWr8f7778Nut+PBBx/Exo0b8fOf/1zT8hgl5JfGdkaMlMCr70RZMwaJCLLNTsh00isgV5jwxFdjtWvXDt9//z3uuecevPrqqxgxYoTf7yXGKCFNtO6Iat+OG4/AMHOMWPWwTV0Mk/oxTIwzb948PPTQQ86XEuuBUUI+0XNHxFkUdVkhSEhdDBNjPP7447ovkye6ktcCvSM6O8jm8kVyskqQyDpLUkO28ckWjdyGqIlRQo0SFQmMFPkwSOQi2zgZJuQvHr6hesn2hObhHrGsEiTkH5nOLyH1cKaEPJItSDzhLIpxrBQkss0+NEa18RqN2we1MErIhao7eR7qCRwGifxkG7dsM1zcJqiDUUJOZnriMlL0wd8dacUwIS0YJWSJHTcjxXeifk+cJdFGxvHLFibFA/jclx1PdLUwK++cedJswxgkapLtTdWIfMUosSgrB4knfKfZ60SuF7L9VU364KtxyBc8fGMxPHzROKse6rFqkJhllqSGjD8Pg5O8xSixCKvtYPVkhUhhkJiLjD8Xw4S8wcM3FmDWHakoZjsfxapBQkTy4UyJiZn5L3uZqDyTYuUgkXE2QU8y/nyiH3OSH6PEpFTbOZqJKoEi+/jIfwwTUg0P35gMdzRykfVQj+j1RPSOScadtZXwFTlUH0aJSYjeyZB3ZIgU0esKg8RYfO8SUgkP35iA6J0MaWf0+Sii1xXRQWJVMoYY1wXyhDMlChO9gyH9BXImRfT6IsNOSMads5XxMA7VxZkSRYnewZAx9JhFkeGkWwaJeLL+/DKsGyQPzpQoRvTOhcTRMosiw/rCnQ4ReYszJYqQ4a9dkktj56PIsL7IEiSyzhIYTdbfgyzrCYnHKFGADDsXkp9sb+Imy45G1h2xKLL+PmRZX0gsRonEZNm5EBEZgWFCjBJJMUZIZbLsXGSdFRCNvxeSFaOEiHQny8s8+aZhnvH3QrLiq2+IKCBqh4nImZO6O2ArzhLIGiGyxCvJg1FCRAFXd+cjS6SYNVAYIaQqRgkRGU6WSDHLLAojhMyCUUJEwvFQj28YIWRWjBIJ8ZU35K/a7/aq2vokyywKIM+hHkYIWQWjhMjk6r4dPSNFG6NmUWQNEIARQoHHKCEymcY+E8fT/6sUKmY71MMIIfoJo4TIRLz5kD5vbqdKpMgyiwJ4HymMEKL6MUqITEJrkHizLEaK72SOD4ABQnJilBBRo1Q95CNTpIjGCCEVMEqITEDPWRKt96lapJg9UBghpCJGCZHiRASJJ6pFitlmURghZAamiJLwDftw9cEBoodBZDhZgsQT1Q75qBYpjBAyI1NECQCEr9uNH3+VInoYRIaROUjqo9JsimyHehghZAWmiRKAYUKkGlUiRcQsCiOErMhUUQIwTESpvZHmxjTwVJwl8YYqh3z0jhQ+Z4iuM12UAAwTo9XdIA+84xA3sgFk1iCpjwqzKb5GCp8fRJ6ZMkoAholR6tv41lzOja++rBYknqgWKYx0Iu8FiR5AIIWv2y16CJYnwwmCZsEg8Sy+wOH2JRMGCZH3TB0lAMMkkLwNDoYJGU3mSCGi+pn28A0Flq+hwcM5/uGO1T8qHPIhIgvMlACcLdGbPzMfnDXxHYNEf5xJIZKTZWZKeOKrPvSICp745z3uMI2hykuRyZWvz49rVx34OkBjIX1YJkoAholMGCaNY5CIxUM++uM6TY2xVJQADBN/6H3oheeZkEqsGCmMCDKa5aIEYJhoEchzQThr4o47A/nJeMiH6w2pzpJRAjBMfGHEyakMk59wx6IuX2dT+FgTubJslAAME28Y+WoZHs7hTsps+HgS+cYSLwluCF8uLB+rvmyYOzAisjrLRwnAMKmPyDiwapgQEVkZo+T/MUxcyRAFA+84JMU4jMBZEiIiRokLhsl1soWAbOPRG4OEiOg6oVGyY8cOjBgxAvHx8bDZbPjwww9FDgcAw0TWAJB1XP5ikBCRnlasWIGkpCSEhYUhJSUFe/bsqfe6r732GgYOHIiWLVuiZcuWGDJkSIPXN4LQKCkvL0fPnj2xYsUKkcNwY/UwkZWVDucQEflq7dq1yMjIQHZ2Ng4cOICePXti6NChOHfunMfrb9++HaNGjcK2bduwa9cuJCQk4Je//CXOnDlj8Mh/IjRKhg0bhueeew6/+tWvRA7DIyuGiSo7fFXG2RjOkhCRnpYuXYrHH38c48ePR9euXZGXl4ef/exneOONNzxe/+9//zueeuop9OrVC126dMHrr78Ou92O/Px8g0f+E6XOKamsrERZWZnLF+lDtR29auOti0FCRN6ou8+rrKz0eL2qqirs378fQ4YMcV4WFBSEIUOGYNeuXV7d15UrV3D16lVER0frMnYtlHrztJycHMyfP1/0MIj8wiAhMr/wDfvQxNZU8+2vOa4CABISElwuz87Oxrx589yuf+HCBVRXVyM2Ntbl8tjYWBw+fNir+5w5cybi4+NdwsZoSkVJZmYmMjIynN+XlZW5PWB6EfVOr6I+O6Pw065KzT5Y+V1ficg6Tp8+jcjISOf3oaGhAbmfRYsWYc2aNdi+fTvCwsICch/eUCpKQkNDA/aA1GbVt55XJUxUDhLOkhAZp+4fefYKG7BB0GA0ioyMdImS+rRu3RrBwcEoKSlxubykpARxcXEN3vbFF1/EokWL8O9//xs9evTwa7z+UipKjCAySER/wqjsVI4RItKO28bGhYSEoHfv3sjPz0daWhoAOE9anTx5cr23e/7557FgwQJs3rwZycnJBo22fkKj5PLlyzh27Jjz+xMnTqCoqAjR0dFITEw0dCxWnR2pS9bZEgYJkTkwMAInIyMD6enpSE5ORt++fZGbm4vy8nKMHz8eADB27Fi0a9cOOTk5AIDFixcjKysLq1evRlJSEoqLiwEAERERiIiIEPIzCI2Sffv24e6773Z+X3O+SHp6Ot58803DxiFDkMj0RJUtTMwSJDx0Q2Yi0zaLrhs5ciTOnz+PrKwsFBcXo1evXti0aZPz5NdTp04hKOinF92uXLkSVVVV+PWvf+2ynPpOpjWC0Ci566674HCI3VDLECRUP7MECZHMGBjmMXny5HoP12zfvt3l+5MnTwZ+QD6y9DklsgSJjBsEGWZLzBQknCUhI8i4LSHyhWWjRJYgkZnIMDFTkBBpwcAgK7JklMgUJNzwuDNbkHCWxNr4HCfynuWiRKYgUYHRsyVmCxKyJoYIkTaWihLZgkSVDZdRYWLGIOEsiXWo8nwmkpklokS2GCF3ZgwSMj+GCJG+TB8lsgaJahuzQM6WMEhIJao9d4lUYuookTVIVBWIMDFzkPDQjXkwRNTQ2Pap6nIVvjZoLKSNaaNE5iDhBu46MwcJqY/PU+OIfk8kkocpo0TmIFGdXrMlZg8SzpKoiSHiPYYEBYLpokT2IDHDRs+fMDF7jJB6zPCc9AVjgmRmqiiRPUiszipBwlkS+akYIowJsgLTRIkKQaLihrA+vs6WWCVISE4in3uMCSLvmSJKfhyRbI4fRDHehomVgoSzJPJgiBCph/tyCigrBQmJxxAhUhujxCBmOnRTW0OzJQwSMgJDhMg8GCXkN09hYsUg4aEb4zBEiMyJUWIAs86S1MeKQUKBxxAhMj9GCemiZrbEqkHCWZLAYIgQWQujJMCsNEti1SAhfTFEqD6rEgv9un3ZpWqs0WksFBiMEiI/cZbEf6JChBHiH38jgaguRkkAWWmWhMhXDBF9MRDIDBglRGQYK4YIY4HIe4ySAOEsiTXw0E3DrHp+CEOESBtGCRHpiiFCRFoxSgKAsyTWwFmSnzBEiEgPjBIi0oQhQkR6Y5TojLMk1mDVWRKGCBEFEqOEiBrEECEiozBKdMRZEmuwwiwJQ4SIRGCU+IkhYj1nB9ksESZEREYLEj0AFZ0dZHN+kTXxsQ8ckZ+hNP7UQGH3TUSMEq/UjhDujKiGWdcFGX4uhgmRNTFK6sEIIW9wHTEnhgmRGIySWhgipBXXGf2JnC0BGCZEIlg+ShgipBczrENm+Bn0xDAhMpYlo4QhQoHCdUpfomdLAIYJkZEsEyUMETIK1zHzYZgQGcPU71PCnQOJouJ7mcj6fCn8tKvQN1SrMf7UQL65miLqi8iqy1UAjhs7GPKJ6aJE1g0rWY+KYUINY5jIg7NX5mSKKCkeYENQGGNElLp/xcpwHoAsGCb6kGW2BGCYGInhYT2miBIyVmM7h4F3HGKY1KJCmHCG0TcME/0wPKg2Rgl5xde/Umuuzzi5ToUwkZ1MsyUAw8QXDA/yFqOE6qXHDoCzJj+RNUw4S6Idw+QnDA/SA6OEnAL1VyhnTX4ia5ioQrbZEsB6YcL4oEBilFickRt4zppcxzAxH7OFCcODRGGUWJDIvzQ5a3KdLGGi4qEbGWdLAPXCxErhUbO9sVdUCB4JNYZRYgEybsA5ayJPmKiIYeIdK4VHbVbftqiMUWJSMm6w62KYMEzMyOgwsWp41GX1bYlZMEpMRIUQqYuHc8SFiYqHbmqTdbYE0D9MGB7urLzNMDNGicJk3SBrYfVZE86YmI+vYcLwaJiVtw9WwihRjJlCpC6rz5rUzFwwTrwn82wJ4B4mDA/vWXU7YHWMEgXIvNENBM6aBH7WRPVDNyphiHjHys95+gmjREJWixBPOGvCwznekn22hDyz6nObGsYokQQ3qp5ZedaEYUJmIvp5HF/gwLWrDnwtdBTUGEaJQAwR71h51iQQYWLGQzecLZGTyOcsg15NjBIDcaPpH6vOmnDGhFQh+vnJ54n6GCUBxhDRl1VnTfQKEzPOktTgbInxRD8PGSHmwygJAG4YA8+KsyacMSHRRD/nuP6bH6NEB4wQMaw4a8IwaRhnS/Ql+rnFdd16GCUaccMnD6vNmmgNEzMfuiF9iH4eMUKIUeIDhoi8rDZrwhmT+nG2xHuiny9ch6kuRkkDuGFTj5VmTRgm9WOYeCb6uSF6fQ3fsE/o/VPjgkQPAABWrFiBpKQkhIWFISUlBXv27BE2loF3HHJ+kZqs9Ph5e0iGh26sq/DTrs4vEeILHM4vo4Wv2+3yZQW+7k/fffdddOnSBWFhYejevTs+/vhjg0bqmfCZkrVr1yIjIwN5eXlISUlBbm4uhg4diiNHjiAmJsaQMVhlB2Y1Vpk14YyJZ1adLRG9zoteF60SH574uj/9z3/+g1GjRiEnJwf33XcfVq9ejbS0NBw4cADdunUT8BMANofDIXQNSklJQZ8+fbB8+XIAgN1uR0JCAqZMmYJZs2Y1eNuysjJERUWhw+LnEBQW5vV9WnFDZWWiN9JGaWhnYNWZEis810Wv3ypFyDXHVWzHRygtLUVkZGRAxlOzX7oLD6CJranm5WgZq6/705EjR6K8vBwbN250XnbHHXegV69eyMvL0zx2fwidKamqqsL+/fuRmZnpvCwoKAhDhgzBrl273K5fWVmJyspK5/elpaUAAHtFhVf3l5p8+Pr9XvZn1KSalG5FAIBd+7qIHUiAfZMCxO1030EUD7AB3j1FTKdg+w3O571ZuK/Hxj64ddexa4beu/t5Ib7c/zVcBQAY8bf4NVwF/LibmrGWlZW5XB4aGorQ0FC36/u6PwWAXbt2ISMjw+WyoUOH4sMPP9Q+cD8JjZILFy6guroasbGxLpfHxsbi8GH3DUlOTg7mz5/vdvnp7Oe8uj9+EBOZncd1fIPRo5ALn/f6MsPv87vvvkNUVFRAlh0SEoK4uDjsLPb/3IyIiAgkJCS4XJadnY158+a5XdfX/SkAFBcXe7x+cXGxfwP3g/BzSnyRmZnpUnUXL15Ehw4dcOrUqYCtYGZWVlaGhIQEnD59OmBTmWbF3512/N35h78/7UpLS5GYmIjo6OiA3UdYWBhOnDiBqqoqv5flcDhgs7keevU0S2ImQqOkdevWCA4ORklJicvlJSUliIuLc7t+fdNWUVFRfHL6ITIykr8/jfi7046/O//w96ddUFBgX3gaFhaGMB/Oc9SDr/tTAIiLi/Pp+kYQ+pLgkJAQ9O7dG/n5+c7L7HY78vPzkZqaKnBkRERE6tCyP01NTXW5PgBs2bJF6P5X+OGbjIwMpKenIzk5GX379kVubi7Ky8sxfvx40UMjIiJSRmP707Fjx6Jdu3bIyckBAEybNg2DBg3CkiVLMHz4cKxZswb79u3Dq6++KuxnEB4lI0eOxPnz55GVlYXi4mL06tULmzZtcjv5xpPQ0FBkZ2eb/hhboPD3px1/d9rxd+cf/v60M/vvrrH96alTp1wOXfXr1w+rV6/GnDlzMHv2bNx000348MMPhb1HCSDB+5QQERERAZK8zTwRERERo4SIiIikwCghIiIiKTBKiIiISApKR4mvH9FM19+qv0+fPmjevDliYmKQlpaGI0eOiB6WkhYtWgSbzYbp06eLHooyzpw5g0cffRStWrVCeHg4unfvjn379jV+Q4urrq7G3Llz0bFjR4SHh+PGG2/Es88+a8hnuKhox44dGDFiBOLj42Gz2dw+y8XhcCArKwtt27ZFeHg4hgwZgi+//FLMYMmFslFS8xHN2dnZOHDgAHr27ImhQ4fi3LlzoocmtYKCAkyaNAmffvoptmzZgqtXr+KXv/wlysvLRQ9NKXv37sUrr7yCHj16iB6KMn744Qf0798fTZs2xSeffIJDhw5hyZIlaNmypeihSW/x4sVYuXIlli9fjv/93//F4sWL8fzzz+Pll18WPTQplZeXo2fPnlixYoXH/3/++eexbNky5OXlYffu3WjWrBmGDh2KCi8/3JUCyKGovn37OiZNmuT8vrq62hEfH+/IyckROCr1nDt3zgHAUVBQIHooyrh06ZLjpptucmzZssUxaNAgx7Rp00QPSQkzZ850DBgwQPQwlDR8+HDHhAkTXC578MEHHaNHjxY0InUAcKxbt875vd1ud8TFxTleeOEF52UXL150hIaGOv7xj38IGCHVpuRMSc1HNA8ZMsR5WWMf0UyelZaWAkBAP6DKbCZNmoThw4e7rH/UuPXr1yM5ORkPPfQQYmJicNttt+G1114TPSwl9OvXD/n5+Th69CgA4PPPP8fOnTsxbNgwwSNTz4kTJ1BcXOzy/I2KikJKSgr3HxIQ/o6uWmj5iGZyZ7fbMX36dPTv31/oO/ipZM2aNThw4AD27t0reijKOX78OFauXImMjAzMnj0be/fuxdSpUxESEoL09HTRw5ParFmzUFZWhi5duiA4OBjV1dVYsGABRo8eLXpoyikuLgYAj/uPmv8jcZSMEtLHpEmT8MUXX2Dnzp2ih6KE06dPY9q0adiyZYvhnwBqBna7HcnJyVi4cCEA4LbbbsMXX3yBvLw8Rkkj3nnnHfz973/H6tWrceutt6KoqAjTp09HfHw8f3dkKkoevtHyEc3kavLkydi4cSO2bduG9u3bix6OEvbv349z587h9ttvR5MmTdCkSRMUFBRg2bJlaNKkCaqrq0UPUWpt27ZF165dXS675ZZbcOrUKUEjUsczzzyDWbNm4eGHH0b37t0xZswYzJgxw/nBauS9mn0E9x9yUjJKtHxEM13ncDgwefJkrFu3Dlu3bkXHjh1FD0kZgwcPxsGDB1FUVOT8Sk5OxujRo1FUVITg4GDRQ5Ra//793V5+fvToUXTo0EHQiNRx5coVlw9SA4Dg4GDY7XZBI1JXx44dERcX57L/KCsrw+7du7n/kICyh28a+4hm8mzSpElYvXo1PvroIzRv3tx5DDUqKgrh4eGCRye35s2bu51706xZM7Rq1Yrn5HhhxowZ6NevHxYuXIjf/OY32LNnD1599VWhH5OuihEjRmDBggVITEzErbfeis8++wxLly7FhAkTRA9NSpcvX8axY8ec3584cQJFRUWIjo5GYmIipk+fjueeew433XQTOnbsiLlz5yI+Ph5paWniBk3XiX75jz9efvllR2JioiMkJMTRt29fx6effip6SNID4PFr1apVooemJL4k2DcbNmxwdOvWzREaGuro0qWL49VXXxU9JCWUlZU5pk2b5khMTHSEhYU5brjhBscf//hHR2VlpeihSWnbtm0et3Pp6ekOh+P6y4Lnzp3riI2NdYSGhjoGDx7sOHLkiNhBk8PhcDhsDgffEpCIiIjEU/KcEiIiIjIfRgkRERFJgVFCREREUmCUEBERkRQYJURERCQFRgkRERFJgVFCREREUmCUEBERkRQYJUQKOHnyJGw2G2w2G3r16iVkDPPmzXOOITc3V8gYiMjcGCVECvn3v//t8kFiRnr66afx7bff8lOliShglP1APiIratWqFVq1aiXkviMiIhAREcFPQyaigOFMCZHBzp8/j7i4OCxcuNB52X/+8x+EhIT4PAsybtw4pKWlYeHChYiNjUWLFi3wpz/9CdeuXcMzzzyD6OhotG/fHqtWrXLepuZQ0DvvvIOBAwciPDwcffr0wdGjR7F3714kJycjIiICw4YNw/nz53X7uYmIGsMoITJYmzZt8MYbb2DevHnYt28fLl26hDFjxmDy5MkYPHiwz8vbunUrzp49ix07dmDp0qXIzs7Gfffdh5YtW2L37t144okn8Lvf/Q7ffPONy+2ys7MxZ84cHDhwAE2aNMEjjzyCP/zhD3jppZdQWFiIY8eOISsrS68fm4ioUYwSIgHuvfdePP744xg9ejSeeOIJNGvWDDk5OZqWFR0djWXLlqFz586YMGECOnfujCtXrmD27Nm46aabkJmZiZCQEOzcudPldk8//TSGDh2KW265BdOmTcP+/fsxd+5c9O/fH7fddhsmTpyIbdu26fHjEhF5heeUEAny4osvolu3bnj33Xexf/9+hIaGalrOrbfeiqCgn/6+iI2NRbdu3ZzfBwcHo1WrVjh37pzL7Xr06OFyGwDo3r27y2V1b0NEFEicKSES5KuvvsLZs2dht9tx8uRJzctp2rSpy/c2m83jZXa7vd7b2Ww2j5fVvQ0RUSBxpoRIgKqqKjz66KMYOXIkOnfujMceewwHDx5ETEyM6KEREQnDmRIiAf74xz+itLQUy5Ytw8yZM3HzzTdjwoQJoodFRCQUo4TIYNu3b0dubi7++te/IjIyEkFBQfjrX/+KwsJCrFy5UvTwiIiEsTkcDofoQRBRw06ePImOHTvis88+E/Y28zWSkpIwffp0TJ8+Xeg4iMh8OFNCpJB+/fqhX79+Qu574cKFiIiIwKlTp4TcPxGZH2dKiBRw7do15yt0QkNDkZCQYPgYvv/+e3z//fcArr8BXFRUlOFjICJzY5QQERGRFHj4hoiIiKTAKCEiIiIpMEqIiIhICowSIiIikgKjhIiIiKTAKCEiIiIpMEqIiIhICowSIiIiksL/AeE3i1SIRvApAAAAAElFTkSuQmCC", "text/plain": [ "

" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "data.plot.contourf()" ] }, { "cell_type": "markdown", "id": "ca0589eb-8d3b-4a1c-bba8-71c528144745", "metadata": {}, "source": [ "## Cirumnavigate return of `xarray.DataArray` objects\n", "\n", "In certain cases, there may be no requirement to return `xarray.DataArray` objects, and it may be more convenient to work with the default interface, hence `numpy.array` objects:" ] }, { "cell_type": "markdown", "id": "560f8b30-bbc4-4f60-b89b-df480508b695", "metadata": {}, "source": [ "If we got the `xarray` object already, just call the property `.values`. Otherwise, we have the following two options to retrieve `numpy.array`:" ] }, { "cell_type": "code", "execution_count": 7, "id": "27f42721-c7d6-40e9-beed-01a918757e2f", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "numpy.ndarray" ] }, "execution_count": 7, "metadata": {}, "output_type": "execute_result" } ], "source": [ "with h5tbx.File(h5.hdf_filename) as h5:\n", " data_np = h5.vel.values[:]\n", "type(data_np)" ] }, { "cell_type": "markdown", "id": "1c3da562-5914-4ed4-9294-d170a5a58a4d", "metadata": {}, "source": [ "Using the configuration setter just for this code snippet (using context manager syntax):" ] }, { "cell_type": "code", "execution_count": 8, "id": "155f9ff4-aa61-4539-882d-d8d3f9348592", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "numpy.ndarray" ] }, "execution_count": 8, "metadata": {}, "output_type": "execute_result" } ], "source": [ "with h5tbx.set_config(return_xarray=False):\n", " with h5tbx.File(h5.hdf_filename) as h5:\n", " data_np = h5.vel.values[:]\n", "type(data_np)" ] }, { "cell_type": "markdown", "id": "99459d8c-44a5-4f3a-b2d2-cf555ccd222a", "metadata": {}, "source": [ "## Selecting data (`.sel`)\n", "\n", "HDF5 datasets may sometimes be very large. Hence it is ineffcient to slice a larger array and then use the useful method of (selecting)[https://docs.xarray.dev/en/stable/user-guide/indexing.html]. The `h5rdmtoolbox` allows to call `.sel` prior to the above slicing, to reduce the data loaded to the RAM:" ] }, { "cell_type": "code", "execution_count": 10, "id": "ef6c4a7e-68d3-49de-909b-b7a57f806dee", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "available coords to select from: dict_keys(['y', 'x'])\n" ] }, { "data": { "text/html": [ "

\n", "

<xarray.DataArray 'vel' (x: 5)>\n",
       "0.3443 0.3647 0.09872 0.8726 0.3015\n",
       "Coordinates:\n",
       "    y        float64 2.0\n",
       "  * x        (x) float64 0.0 2.5 5.0 7.5 10.0\n",
       "Attributes:\n",
       "    long_name:  velocity\n",
       "    units:      m/s

\n", "\n", "

\n", " \n", " \n", " /(4)\n", " \n", "
\n", "\n", "
\n", "\n", "
\n", "\n", "
\n", "\n", "
\n", "

\n", "

" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "h5tbx.dump(h5)" ] }, { "cell_type": "markdown", "id": "ccc401dc-9dde-4f5f-b04a-26b681fbfed5", "metadata": {}, "source": [ "The ancillary dataset will appear as a `xarray` coordinate when the dataset is sliced:" ] }, { "cell_type": "code", "execution_count": 14, "id": "6de1dffc-f602-4e42-b7d9-0993ce3e3c6f", "metadata": {}, "outputs": [], "source": [ "with h5tbx.File(h5.hdf_filename) as h5:\n", " u = h5.vel[()]" ] }, { "cell_type": "code", "execution_count": 15, "id": "2ce04d0c-b674-4119-baff-917153645778", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "available ancillary datasets: {'uncertainty': }\n" ] }, { "data": { "text/plain": [ "Coordinates:\n", " y float64 3.0\n", " * x (x) float64 0.0 2.5 5.0 7.5 10.0\n", " uncertainty (x) float64 0.02477 0.02484 0.02294 0.02668 0.02445" ] }, "execution_count": 15, "metadata": {}, "output_type": "execute_result" } ], "source": [ "with h5tbx.File(h5.hdf_filename) as h5:\n", " print('available ancillary datasets: ', h5.vel.ancillary_datasets)\n", " data = h5.vel.sel(y=3.1, method='nearest')\n", "data.coords" ] }, { "cell_type": "markdown", "id": "78cdd07c-6d08-4d83-abd2-d26ba95c39a2", "metadata": {}, "source": [ "## Unit interface\n", "\n", "It is reasonable to provide physical units to the datasets. This is probably one of the most important meta information about a dataset. You might have noticed its usage in the examples above.\n", "\n", "If the dataset is sliced as is, the returned `xarray.DataArray` will have the respective unit. In some cases, it might be necessary to convert the data (and/or its coordinates) into another unit. The `h5rdmtoolbox` provides an interface for this purpose. Address the desired dataset and call `to_units(...)` on it, which will call the interface class. The latter implements all major methods, like `sel`, `isel` and `__getitem__`.\n", "\n", "Note, that this feature must be \"enabled\" by importing the `units` modue from `extensions`. The implementation principle is discussed [here](../misc/Extensions.ipynb)." ] }, { "cell_type": "code", "execution_count": 16, "id": "42e207df-426a-431d-ab1b-de953eda1144", "metadata": {}, "outputs": [], "source": [ "from h5rdmtoolbox.extensions import units" ] }, { "cell_type": "code", "execution_count": 17, "id": "8535f777-4725-4d0b-b0f5-1dcca5d1026b", "metadata": {}, "outputs": [ { "data": { "text/html": [ "

\n", "

<xarray.DataArray 'vel' (y: 11, x: 5)>\n",
       "583.2 734.7 389.9 701.3 79.84 317.7 ... 484.6 546.2 70.77 386.4 299.8 247.3\n",
       "Coordinates:\n",
       "  * y            (y) float64 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0\n",
       "  * x            (x) float64 0.0 0.0025 0.005 0.0075 0.01\n",
       "    uncertainty  (y, x) float64 0.02454 0.02421 0.02566 ... 0.02507 0.0234\n",
       "Attributes:\n",
       "    ANCILLARY_DATASETS:  ['uncertainty']\n",
       "    long_name:           velocity\n",
       "    units:               mm/s

" ], "text/plain": [ "\n", "583.2 734.7 389.9 701.3 79.84 317.7 ... 484.6 546.2 70.77 386.4 299.8 247.3\n", "Coordinates:\n", " * y (y) float64 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0\n", " * x (x) float64 0.0 0.0025 0.005 0.0075 0.01\n", " uncertainty (y, x) float64 0.02454 0.02421 0.02566 ... 0.02507 0.0234\n", "Attributes:\n", " ANCILLARY_DATASETS: ['uncertainty']\n", " long_name: velocity\n", " units: mm/s" ] }, "execution_count": 17, "metadata": {}, "output_type": "execute_result" } ], "source": [ "with h5tbx.File(h5.hdf_filename) as h5:\n", " data = h5.vel.to_units('mm/s', x='m')[()]\n", "data" ] }, { "cell_type": "markdown", "id": "06379937-9d96-40b1-8c3b-df2dfad6f2fb", "metadata": {}, "source": [ "As a side note, `xarray` works well with `pint`, which is a units-library. Please refer to the [pint-xarray documentation](https://#pint-xarray.readthedocs.io/en/stable/) or [to this blog entry about it](https://xarray.dev/blog/introducing-pint-xarray) for detailed explanation and #examples on how it is done.\n", "\n", "The following shows an example usage of it (which takes place under the hood of `to_units()` implemented in the `h5rdmtoolbox`:" ] }, { "cell_type": "code", "execution_count": 18, "id": "dff826e5-35e4-45bc-8bbd-490e15d6deff", "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'mm/s'" ] }, "execution_count": 18, "metadata": {}, "output_type": "execute_result" } ], "source": [ "# note, that we need to call the \"pint\"-accessory every time to make the xarray object\n", "quantified_data = data.pint.quantify(unit_registry=h5tbx.get_ureg())\n", "\n", "u_mm_s = quantified_data.pint.to('mm/s').pint.dequantify()\n", "u_mm_s.units" ] }, { "cell_type": "code", "execution_count": null, "id": "b302e1b4-5efa-4126-bdcc-9b551fd99564", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.19" } }, "nbformat": 4, "nbformat_minor": 5 }