{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dfsu2D\n",
       "Number of elements: 884\n",
       "Number of nodes: 529\n",
       "Projection: UTM-29\n",
       "Items:\n",
       "  0:  Surface elevation <Surface Elevation> (meter)\n",
       "  1:  U velocity <u velocity component> (meter per sec)\n",
       "  2:  V velocity <v velocity component> (meter per sec)\n",
       "  3:  Current speed <Current Speed> (meter per sec)\n",
       "Time: 9 steps with dt=9000.0s\n",
       "      1985-08-06 07:00:00 -- 1985-08-07 03:00:00"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from mikeio import Dfsu\n",
    "\n",
    "filename = \"../tests/testdata/HD2D.dfsu\"\n",
    "dfs = Dfsu(filename)\n",
    "\n",
    "dfs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<mikeio.Dataset>\n",
       "Dimensions: (9, 884)\n",
       "Time: 1985-08-06 07:00:00 - 1985-08-07 03:00:00\n",
       "Items:\n",
       "  0:  Surface elevation <Surface Elevation> (meter)\n",
       "  1:  Current speed <Current Speed> (meter per sec)"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds = dfs.read([\"Surface elevation\",\"Current speed\"]) # to read some variables\n",
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████████████████████| 9/9 [00:00<00:00, 14745.60it/s]\n"
     ]
    }
   ],
   "source": [
    "Dfsu.show_progress = True # Turn on progress bar, useful to get progress for long-running tasks\n",
    "ds = dfs.read()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Surface elevation</th>\n",
       "      <th>U velocity</th>\n",
       "      <th>V velocity</th>\n",
       "      <th>Current speed</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>count</th>\n",
       "      <td>7956.000000</td>\n",
       "      <td>7956.000000</td>\n",
       "      <td>7956.000000</td>\n",
       "      <td>7956.000000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>mean</th>\n",
       "      <td>0.034152</td>\n",
       "      <td>-0.027504</td>\n",
       "      <td>0.051639</td>\n",
       "      <td>0.073658</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>std</th>\n",
       "      <td>0.534890</td>\n",
       "      <td>0.088108</td>\n",
       "      <td>0.092440</td>\n",
       "      <td>0.119606</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>min</th>\n",
       "      <td>-0.758378</td>\n",
       "      <td>-1.280627</td>\n",
       "      <td>-0.811804</td>\n",
       "      <td>0.000082</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>25%</th>\n",
       "      <td>-0.579662</td>\n",
       "      <td>-0.022283</td>\n",
       "      <td>0.019642</td>\n",
       "      <td>0.029917</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>50%</th>\n",
       "      <td>0.101636</td>\n",
       "      <td>-0.011019</td>\n",
       "      <td>0.041748</td>\n",
       "      <td>0.046848</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>75%</th>\n",
       "      <td>0.462013</td>\n",
       "      <td>0.004649</td>\n",
       "      <td>0.062974</td>\n",
       "      <td>0.070259</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>max</th>\n",
       "      <td>1.172707</td>\n",
       "      <td>0.056265</td>\n",
       "      <td>1.518996</td>\n",
       "      <td>1.714988</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "       Surface elevation   U velocity   V velocity  Current speed\n",
       "count        7956.000000  7956.000000  7956.000000    7956.000000\n",
       "mean            0.034152    -0.027504     0.051639       0.073658\n",
       "std             0.534890     0.088108     0.092440       0.119606\n",
       "min            -0.758378    -1.280627    -0.811804       0.000082\n",
       "25%            -0.579662    -0.022283     0.019642       0.029917\n",
       "50%             0.101636    -0.011019     0.041748       0.046848\n",
       "75%             0.462013     0.004649     0.062974       0.070259\n",
       "max             1.172707     0.056265     1.518996       1.714988"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds.describe()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Find which element is nearest to POI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "idx = dfs.find_nearest_elements(606200, 6905480)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract a subset of the dataset from this element. (Discrete values, no interpolation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<mikeio.Dataset>\n",
       "Dimensions: (9,)\n",
       "Time: 1985-08-06 07:00:00 - 1985-08-07 03:00:00\n",
       "Items:\n",
       "  0:  Surface elevation <Surface Elevation> (meter)\n",
       "  1:  U velocity <u velocity component> (meter per sec)\n",
       "  2:  V velocity <v velocity component> (meter per sec)\n",
       "  3:  Current speed <Current Speed> (meter per sec)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "selds = ds.isel(idx=idx)\n",
    "selds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert to a dataframe, for convenience."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>Surface elevation</th>\n",
       "      <th>U velocity</th>\n",
       "      <th>V velocity</th>\n",
       "      <th>Current speed</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>1985-08-06 07:00:00</th>\n",
       "      <td>0.459460</td>\n",
       "      <td>0.006372</td>\n",
       "      <td>-0.007143</td>\n",
       "      <td>0.009572</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1985-08-06 09:30:00</th>\n",
       "      <td>0.806965</td>\n",
       "      <td>0.010517</td>\n",
       "      <td>0.003438</td>\n",
       "      <td>0.011064</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1985-08-06 12:00:00</th>\n",
       "      <td>0.100285</td>\n",
       "      <td>0.011300</td>\n",
       "      <td>0.012926</td>\n",
       "      <td>0.017169</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1985-08-06 14:30:00</th>\n",
       "      <td>-0.727009</td>\n",
       "      <td>0.010402</td>\n",
       "      <td>0.010033</td>\n",
       "      <td>0.014452</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1985-08-06 17:00:00</th>\n",
       "      <td>-0.579541</td>\n",
       "      <td>0.007293</td>\n",
       "      <td>-0.000902</td>\n",
       "      <td>0.007349</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "                     Surface elevation  U velocity  V velocity  Current speed\n",
       "1985-08-06 07:00:00           0.459460    0.006372   -0.007143       0.009572\n",
       "1985-08-06 09:30:00           0.806965    0.010517    0.003438       0.011064\n",
       "1985-08-06 12:00:00           0.100285    0.011300    0.012926       0.017169\n",
       "1985-08-06 14:30:00          -0.727009    0.010402    0.010033       0.014452\n",
       "1985-08-06 17:00:00          -0.579541    0.007293   -0.000902       0.007349"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "df = selds.to_dataframe()\n",
    "df.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<AxesSubplot:>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEECAYAAAA2xHO4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABN90lEQVR4nO3deXxU1fn48c8zSxJIQtgChCQQdgSyQBYQZV9VFEFWrYIbUnfbWvFrf2qttlq1WltbZRMXlE3EtYqACAgCCQbZ90ACYQkQSAJJJjPn98dMQghJCGQmdyY579drXnPnzrn3PjeZO8+ce+85R5RSaJqmaVpFTEYHoGmapnk3nSg0TdO0SulEoWmaplVKJwpN0zStUjpRaJqmaZXSiULTNE2rlMXoACrTtGlTFRUVZXQYmqZpPiMlJSVLKRXqznV6daKIiooiOTnZ6DA0TdN8hogcdPc69aknTdM0rVI6UWiapmmV0olC0zRNq5ROFJqmaVqldKLQNE3TKuWWRCEiw0Vkl4jsFZFp5bwfIiJfishmEdkmIne7Y7uapmma51U7UYiIGXgbuAHoAkwUkS5lij0EbFdKxQL9gddFxK+62/YGhUUOo0PQNE3zKHfUKJKAvUqp/UqpQmAeMLJMGQUEi4gAQcApoMgN2zbUOz/uo/sLS9mRedboUDRN0zzGHYkiHEgv9TrDNa+0fwPXAEeALcBjSimf/in+0c8Hefl/O8krtPPuj/uMDkfTNM1j3JEopJx5ZYfNGwakAi2BOODfItKg3JWJTBGRZBFJPnHihBvCc7/PUw/z/z7fyqDOzZjcO4ovf80k4/Q5o8PSNE3zCHckigwgstTrCJw1h9LuBhYrp73AAaBzeStTSk1XSiUopRJCQ93aXYlbrNh5jN8v2ExSVGPevqMHU/q2RYDZa9KMDk3TNM0j3JEoNgIdRKSN6wL1BOCLMmUOAYMARKQ50AnY74Zt16if95/ktx9tokvLBsyclECA1UzLhvW4Ja4l8zYeIvtcodEhapqmuV21E4VSqgh4GPgO2AEsUEptE5GpIjLVVewvQG8R2QIsB55SSmVVd9s1aUvGGe57P5nIxvWZc3cSwQHWkvem9G3LuUI7H/3s9r64NE3TDOeW3mOVUt8A35SZ906p6SPAUHdsywh7j+cw6b0NhNSz8uG9STQOvPjO3s4tGtC/Uyhz1qZxX5+2BFjNBkWqaZrmfrpl9mWknzrHb2ZuwCTC3Pt6EhZSr9xyD/RtR1ZuIYs3Ha7hCDVN0zxLJ4pKHM/J585Z6zlXWMSH9yYR1TSwwrK92jYmJiKEGav3Y3eUvelL0zTNd+lEUYEz52zcNWsDx84W8N7dSVwTVu7dvCVEhAf6tuNAVh7fbz9WQ1FqmqZ5nk4U5ThXWMTdczaw/0Qe0++KJ751oyotN7xbC1o1rs87P+5DKV2r0DStdtCJooyCIjsPfJhCano2b02Mo0+HqrflMJuE+/u0ITU9m41ppz0YpaZpWs3RiaKUIruDx+elsnpPFq/cFsPwbmFXvI4x8ZE0DvRj+irdrYemabWDThQuSin+77Mt/G/rUf7fiC6MTYi8/ELlqOdn5q5rW7Nsx3H2HMtxc5Sapmk1TycKnEnixa93sCA5g0cHdeDe69tUa313XRtFgNXE9FU+1/hc0zTtEjpRAP9asZdZaw4wuXcUTwzuUO31NQ70Y3xCJEtSD3PsbL4bItQ0TTNOnU8Uc346wD++381tPSJ4dkQXnENmVN99fdpidyhm/3TALevTNE0zSp1OFIs3ZfD8l9sZ2qU5r9wWjcnkniQBENm4PjdGh/Hxz4fIybe5bb2apmk1rc4miqXbjvLkol+5rn0T3prYHYvZ/X+KB/q2I6egiE82HHL7ujVN02pKnUwUa/dm8fDHvxAdHsL0OxM81olfdEQIvds1YfaaND22tqZpPqvOJYrU9Gzu+yCZNk0DmXN3IoH+bulAt0IP9GvH0bP5fLG57FhOmqZpvqFOJYpdR3OY/N4Gmgb58+G9STSs73f5haqpb4emdG4RzPRV+3DozgI1TfNBdSZRHDp5jjtnrcffYmLufT1p1iCgRrYrIjzQry27j+WycvfxGtmmpmmaO9WJRHHsbD53zPqZQruDD+/tSWTj+jW6/RExLWkZEsC7P+oGeJqm+R63JAoRGS4iu0Rkr4hMq6BMfxFJFZFtIvKjO7ZbFafzCrlz1npO5Rby/t1JdGweXFObLmE1m7i3T1vWHzjFL4d0Z4GapvmWaicKETEDbwM3AF2AiSLSpUyZhsB/gFuUUl2BsdXdblXkFhQxec5G0k6eY8akBGIjG9bEZss1ITGSBgEW3a2Hpmk+xx01iiRgr1Jqv1KqEJgHjCxT5nZgsVLqEIBSyuMn6/Ntdu5/P5mth8/w9u096N2uqac3WalAfwt3Xtuab7cdJS0rz9BYNM1b7Dqaw/A3V7F021GjQ9Eq4Y5EEQ6kl3qd4ZpXWkegkYisFJEUEbmropWJyBQRSRaR5BMnTlxVQEV2B4988gvr9p/ktbExDOnS/KrW426TekdhNZuYsVrXKjRtz7Ecbp/xMzuP5vCXr7frtkZezB2Jorx+L8reB2oB4oGbgGHA/xORjuWtTCk1XSmVoJRKCA2t+qBBxRwOxR8X/cr324/xwsiujOoeccXr8JRmwQHc1iOchSkZZOUWGB2Ophlm7/FcJs5Yj8kkPHdzF9JPnWd+cvrlF9QM4Y5EkQGUHrwhAijbuiwD+FYplaeUygJWAbFu2PZFlFK88NV2Fv9ymN8P6chd10a5exPVdl+fttjsDj5Ym2Z0KJpmiANZedw+42dA8cn9PZncO4rEqEb8a/ke8m12o8PTyuGORLER6CAibUTED5gAfFGmzOdAHxGxiEh9oCewww3bvsgb3+9mzto07u/ThocHtnf36t2iXWgQQ7s05/11B8krKDI6HE2rUQdP5jFx+s8UORQf39+L9s2CERGeHNaZ4zkFfLAuzegQtXJUO1EopYqAh4HvcH75L1BKbRORqSIy1VVmB/At8CuwAZiplNpa3W2XNnP1ft5asZfxCZH8343XuK27cE+Y0rcdZ87bWKCr2lodkn7qHBOn/0xBkZ259/W86Fb1pDaN6dcxlP+s3Kd7W/ZCbmlHoZT6RinVUSnVTin1kmveO0qpd0qVeVUp1UUp1U0p9aY7tltswcZ0Xvx6BzdGt+Cvo6O9OkkAxLduRGJUI2auPkCRXV/A02q/jNPnmDD9Z/IK7Xx0X0+uCWtwSZk/DO1E9jkbM1frMVy8jc+3zP7flkymLf6VPh2a8sb4OMxuHFPCk6b0bcfh7PN8vSXT6FA0zaOOZJ9n4oyfycm3Mfe+nnRtGVJuueiIEG6MbsHM1fs5lVdYw1FqlfHpRLFq9wkenfcL3Vs14t074/G3eKa7cE8Y1LkZ7UIDeffH/SilOwvUaqfMM84kkX3Oxof39qRbePlJotjvhnTkvM3Of1furaEItarw2USRcvAUD3yYQvtmwcyenEh9P892F+5uJpPwQN92bM88y097Txodjqa53bGz+dw+Yz0ncwv54J6kKvWM0L5ZMKN7RPD+uoNknjnv+SC1KvHJRLEj8yx3v7eRFiEBfHBPEiH1rEaHdFVGdm9Js2B/3l21z+hQNM2tjufkM3HGzxw/m8/79yTSvVWjKi/72KAOKKX41wpdq/AWPpcoDmTlceesDQT6W/jw3iRCg/2NDumq+VvM3H1dG1bvyWLr4TNGh6NpbpGVW8DtM9Zz9Ew+c+5JIr514ytaPrJxfSYmtWLBxnQOntTd3XgDn0oUmWfO85uZ63EoxYf39iSiUc12F+4Jt/dsRZC/RXfrodUKJ3MLuGPGejJOn2P25EQSo64sSRR7eEB7LGbhzWV73ByhdjV8JlGczC3gNzPXc/a8jQ/uSaJ9syCjQ3KLkHpWJiZF8tWvmaSfOmd0OJp21U7nFXLHzPWkncxj9qREerVtctXratYggMm927Ak9TC7jua4MUrtavhEojibb2PSexvIOH2emZMSLnvnhK+55/o2CDBrjb5/XPNN2eecSWJ/Vh4zJyXQu331e2ue2q8tQX4WXl+6yw0RatXh9Yki32bnvveT2ZmZwzu/iadnNX6leKuwkHqMjAtn/sZ0Tuv7xzUfc+a8jTtnbWDv8Vym3xlPnw5X3plneRrW92NK37Ys3X6M1PRst6xTuzpenSgU8ODcTWxMO8Ub4+MY0LmZ0SF5zJS+bTlvs/PRzweNDkXTquxsvo27Zq1n59GzvHNnD/p3cu8xevf1bWgS6Mdr3+lahZG8OlFknDrHip3HeenWaG6ObWl0OB7VqUUwAzqF8v66NN2DpuYTcvJtTJq9gW1HzvKfO+IZ2Nn9474E+Vt4cEB71uzNYu3eLLevX6sar04U2edtPDW8M7f3bGV0KDXigX7tyMot5NNNGUaHommVyi0oYvJ7G9mScYZ/397Do4OD3dGzFWEhAby6dJfuxcAgXp0oQoP9+W3/dkaHUWN6tmlMbEQIM1btx+7QB4Tmnc4VFnHPextJTc/mrYndGd6thUe3F2A18+igDvxyKJsVOz0+irJWDq9OFC0aBBgdQo0SER7o1460k+f4frseQ1jzPucL7dwzZyPJB0/x5vg4bowOq5HtjomPIKpJfV79bhcO/SOqxnl1oqiLhnVtQesm9fmv7ixQ8zL5Njv3fbCRDQecN5fU5HVDq9nEE0M6svNoDl/pHpdrnE4UXsZsEu7r05bN6dlsOHDK6HA0DXAmifs/SGbtvpO8OiaWkXHhNR7DzTEt6dwimH8s3YVNj+NSo3Si8EJj4yNoEujH9FW6Ww/NeAVFdqZ+lMLqPVm8clsMt8VHGBKHyST8YWgn0k6e49MUfcNHTXJLohCR4SKyS0T2isi0SsoliohdRMa4Y7u1VYDVzF3XRrF853F2H9PdF2jGKSxy8OBHm1i56wR/Gx3NuIRIQ+MZdE0zurdqyD+X79G3kdegaicKETEDbwM3AF2AiSLSpYJyr+AcW1u7jLuubU09q1nXKjTDFBY5eOjjTSzfeZwXb+3GxCTjb1MXEZ4c1onMM/nMXX/I6HDqDHfUKJKAvUqp/UqpQmAeMLKcco8AnwL6/rYqaBTox/jESD5PPczRM/lGh6PVMTa7g0c/+YXvtx/jz7d05Te9WhsdUone7ZpyXfsm/OeHveQVFBkdTp3gjkQRDqSXep3hmldCRMKBUcA7l1uZiEwRkWQRST5x4oQbwvNd917fBrtD8d5PurNAreYU2R08Pi+Vb7cd5dkRXZjUO8rokC7xh6GdOJlXqI+NGuKORCHlzCt7X+ebwFNKqcueVFRKTVdKJSilEkJD3dO5mK+KbFyfm2JaMnf9Ic7m24wOR6sDiuwOnliwma+3ZPKnm67hnuvbGB1Subq3asSQLs15d9V+ss/pjjQ9zR2JIgMofYUrAjhSpkwCME9E0oAxwH9E5FY3bLvWe6BvW3ILivhEn4/VPMzuUDy56Fe+3HyEaTd05r4+bY0OqVK/H9qR3IIi3tXX8TzOHYliI9BBRNqIiB8wAfiidAGlVBulVJRSKgpYBDyolFrihm3Xet3CQ7iufRNm/3SAgiJ9l4fmGQ6H4o+LfuWzXw7z5LBOTO3n/V3ndG7RgJGxLXnvpwMcz9HX8Typ2olCKVUEPIzzbqYdwAKl1DYRmSoiU6u7fg0e6NuOY2cL+Dy1bEVN06rP4VA8vXgLn27K4InBHXloQHujQ6qyxwd3pMiueHvFXqNDqdXc0o5CKfWNUqqjUqqdUuol17x3lFKXXLxWSk1WSi1yx3brij4dmnJNWANmrNqv+7nR3MrhUDyzZCvzk9N5dGB7HhvcweiQrkhU00DGJUby8YZDeihhD9Its32AiPBA37bsOZ7LD7v03cWaeyilePaLrXyy4RAP9m/HE0M6Gh3SVXl0YAdEhH8u32N0KLWWThQ+4qaYMMIb1uPdH/WFO636lFL8+cvtfPTzIR7o25Ynh3VCpLwbGL1fi5AA7urVmsWbMth7XPdk4Ak6UfgIq9nEvde3YUPaKTYdOm10OJoPU0rxl692MGdtGvde34ZpN3T22SRR7Lf921HPauaN73WtwhN0ovAh4xMjCalnZbquVWhXSSnF3/63k9k/HWBy7yj+dNM1Pp8kAJoE+XNvn7Z8vSWTrYfPGB1OraMThQ8J9LdwZ6/WfLf9KPtP5BodjuZjlFL8/btdTF+1nzt7tea5m7vUiiRR7L4+bWhY38prS3cZHUqtoxOFj5nUOwqr2cTMNbrrAu3KvPH9bv67ch+392zFn2/pWquSBECDACu/7deOlbtO6LFc3EwnCh8TGuzPbT0iWJSSwYmcAqPD0XzEP5ft4a0VexmfEMmLI7thMtWuJFHsrmujaBbsz6vf7dQjRLqRThQ+6P4+bbDZHXywLs3oUDQf8O8Ve3hj2W7GxEfwt9HRtTZJANTzM/PIoA5sTDvNj7vrdqei7qQThQ9qGxrE0C7N+WDdQd3Nslap/67cx2tLdzOqeziv3BZTq5NEsfEJkUQ2rser3+3SDVTdRCcKH/VAv3acOW9j/sb0yxfW6qSZq/fzyrc7uSW2Ja+NjcVcB5IEgJ/FxOODOrLtyFm+3XbU6HBqBZ0ofFSPVo1IimrMrDUH9EDz2iVO5hbwyrc7GdKlOf8YV3eSRLFbu4fTvlkQry/dhV3XKqpNJwofNqVvWw5nn+ebLZlGh6J5mSWpR7DZFU8O64TFXPcOc7NJ+MPQjuw7kcdnvxw2OhyfV/c+QbXIwM7NaN8siHd+3K/v8NBKKKVYmJxObEQIHZsHGx2OYYZ1bUF0eAhvfL9bd9FfTTpR+DCTSZjSty07Ms+yZm+W0eFoXmLr4bPsPJrD2ITIyxeuxUSEJ4d14nD2eX0tr5p0ovBxI+Na0izYX3cWqJVYmJKOv8XEzbEtjQ7FcH06NKVnm8a8tXwv5wr1HYJXSycKH+dvMXPP9W1YszdL93GjkW+zs+SXwwzv1oKQelajwzFcca0iK7eA99ceNDocn+WWRCEiw0Vkl4jsFZFp5bx/h4j86nqsFZFYd2xXc7q9ZyuC/C1M12MH13nfbz/G2fwixsbX7dNOpSVENWZAp1De+XEfZ87bjA7HJ1U7UYiIGXgbuAHoAkwUkS5lih0A+imlYoC/ANOru13tggYBVm7v2Yqvt2TqUb7quAXJ6YQ3rEfvdk2MDsWr/H5oJ86ctzFztf4xdTXcUaNIAvYqpfYrpQqBecDI0gWUUmuVUsWDKPwMRLhhu1opd18XhUlglu4ssM46kn2eNXuzuC0+ok60wL4S3cJDuCkmjFlrDpCVq/tIu1LuSBThQOlbCjJc8ypyL/A/N2xXKyUspB4j48KZvzGd03mFRoejGWDxpgyUgrHx+ndYeX43pCP5Njv/XbnP6FB8jjsSRXk/Xcq9qV9EBuBMFE9VuDKRKSKSLCLJJ07oTr2uxJS+bTlvs/Phz/qiXV2jlGJhSga92jYmsnF9o8PxSu1CgxgTH8GHPx/kSPZ5o8PxKe5IFBlA6StnEcCRsoVEJAaYCYxUSp2saGVKqelKqQSlVEJoaKgbwqs7OjYPZmDnZry/No18m25gVJdsOHCKgyfPMa6Ot524nEcHdQAF/1qhh0y9Eu5IFBuBDiLSRkT8gAnAF6ULiEgrYDFwp1Jqtxu2qVXggb5tOZlXyKKUDKND0WrQwpQMgvwt3NAtzOhQvFpEo/rc3rMVC5IzOJCVZ3Q4PqPaiUIpVQQ8DHwH7AAWKKW2ichUEZnqKvYs0AT4j4ikikhydberlS+pTWNiIxsyY/V+3RlaHZFbUMTXv2Zyc2wY9fzMRofj9R4a0B4/s4k3vte/WavKLe0olFLfKKU6KqXaKaVecs17Ryn1jmv6PqVUI6VUnOuR4I7tapcSEab2bcvBk+f4TnexXCd882sm5212xui2E1USGuzP3ddF8cXmI+zIPGt0OD5Bt8yuhYZ2bUFUk/q8++M+3VlgHbAgOZ22oYH0aNXQ6FB8xgN92xEcYOH1pbpWURU6UdRCZpNwX5+2bM44w3o9yHyttv9ELskHTzMuIRIR3XaiqkLqW5narx3Ldhxj06HTl1+gjtOJopYaEx9Bk0A/3v1R3zNemy1KycBsEkZ3r6zpklaeyb2jaBrkx2vf7TI6FK+nE0UtFWA1M6l3FD/sOsGuozlGh6N5gN2h+HRTBv07htKsQYDR4ficQH8LDw1oz9p9J/lJd9NfKZ0oarE7e7WmntWsOwuspVbtOcGxswWMTdAtsa/W7T1b0TIkgL9/t0tfz6uEThS1WKNAP8YnRvLF5sNkntEtUWubhcnpNA70Y2Dn5kaH4rP8LWYeH9yRzenZfL/9mNHheC2dKGq5e69vg0PBez+lGR2K5kan8wpZtv04t8aF42fRh3F1jO4RTtumgby+dLdue1QB/Qmr5SIb1+em6DA+Xn+I3AI9wldt8XnqYQrtDn3ayQ0sZhNPDOnIrmM5fLn5kt6HNHSiqBMm9Y4it6CIr/RBUGssSM6gW3gDrglrYHQotcJN0WFcE9aAN5btxmZ3GB2O19GJog7o0aohHZsHMU8PMF8rbDtyhu2ZZ3UHgG5kMglPDuvIwZPnWJis+0krSyeKOkBEGJ/YitT0bHYe1V0W+LqFyRn4mU3cEtvS6FBqlQGdmhHfuhFvLd+je18uQyeKOmJ093D8zCbmbdC1Cl9WUGRnSephhnZtTsP6fkaHU6uICE8O68TRs/l8pMd0uYhOFHVEo0A/hnVrwWe/HNa/lnzY8h3HyT5nY6w+7eQRvdo2oU+Hprz9w15y8m1Gh+M1dKKoQyYmRnLmvE33KuvDFiSnExYSwPXtmxodSq315LBOnD5nY/aaNKND8Ro6UdQhvdo2oVXj+nyy4ZDRoWhX4eiZfFbtPsFtPSIwm3QHgJ4SE9GQ4V1bMGP1fj3+vItOFHWIySSMT4zk5/2n9OhePmjxLxk4lLPDR82zfje0I3mFRbyjO9UEdKKoc8bGO3+Nzte3yvoUpRQLkzNIimpMVNNAo8Op9To2D2ZUXDjvr0vj2Nl8o8MxnMUdKxGR4cA/ATMwUyn1cpn3xfX+jcA5YLJSapM7tl3bKaWgqAhV6nHRa1sRqsgGdvvFr0vK2J2vHQ4QEwEm4X5TOns/3012wFEsFjOICTEJmEzOabPJNS2IyXRhvknAbC4z3zVtNjvHQyg732S6UNZkuqiMiGt9xdNahVIOnuZAVh4P9m9ndCg+STkczmPHZkOVOlawV3Rs2XmkSS4HM3fx6X+zuDMxAmUvuvD5LXNsIDhfiwlMZY4PuXBsVfZeyTFYfDxc7XseUO1EISJm4G1gCJABbBSRL5RS20sVuwHo4Hr0BP7req6UPTub7CVLQAFKOb/slB0cdpTd+YzDDsoBDgfKUeSa53DNc04rZQd76fmueQ7HReVVqXVdmHdhmdLvKVVqXXbnA/uF6XJfOxwou3J+eTsUqsiBw+6avmQ5BQ6FeKDvmZGu58yVbl91tSgBZRIwm1AWE5gFsZjBbEIsZkxmM2KxYLJaMJnNmK1WTBYzYrUiFgticb4vVkvJMxYLYrE651mtrrJWxM+1jNUKFj/na6ufq0ypZ38/xOKH+PmD1Yr4+zvnm51jU5d0OFo8cWFGqddSZpYqfxmRUoUuXu/33+2kk+0EQ4NyKNy/y7nOS7Z1aSxKOcBe5PzsFbmeHXYoKnIdR0Wuz5vrtcN1rNjtzs++vdS84mPMfuG5ZH2uz3HxMen8vF94dpZzXHguvazD4fzyttuhyI4qcsbjKLLjKHl2frk77M51qiK7a5mLjzXnehVid4Bdgd2BOBRylYfRiwBrIOO9q1u+tnBHjSIJ2KuU2g8gIvNwfheVThQjgQ+U8wj5WUQaikiYUiqzshXbMg6TOe1pN4R49RSq5DgH55cZ4jomXc8Ok+thBrvp4keR+cKzzSQUmcBmAZsJisxyoawZispMO0xQZBLsrnJFFaxfTK4fGKKcP1pMyvkDw/WD3SzOaaVcOdcBDlcCEtdOlHznOLj4dZlpKX5d/MA5z/laLppfXN75LAhgcs0zOVzPrteiFCbXa4vd4Xw4wGy3YbGD2QEWu/NhtoOlECx5CosdrK75VjtYHKqknMlx4dmkfLvGcqvrceRrY+NwB4colIBDXB8ZcR4X9uJjpfRn/KJjobiM67ixXFimeLmLyxU/xPkosy5lUs6T78UfTNePfgTEDCLOeUoUZhFMZlBIyee/+HgqPV36NaXfdy0n6tIycOH5QhlxHRdlj5NS05T/Ph4Yh8kdiSIcKH3CO4NLawvllQkHKk0UjhaB5ExLogBFoQiFAoU4KDCBDUU+ChsOCnBQIM7nQuWgADuFKPJVEYUUka8cFKgiClQRhcrOeWWjUNmwKYfrg+rcniqTALhMNc7P5Ie/2Q8/s7/z2eSHn7n44e96bcXf7I/V7Ie/2f/C+yY/13wr9UvKOh8VlfMvp5zVZL2q6ubrS3fx9g97WfPUQFo2rHfFy1+WuvhIUA4HDmXHgfMXpV05f6k6KJ524FAObPYCCu35FBQ5H/mu5wK761E8315Agb2QAnsBBfYC8u35FNoLybcXUuhwPhfYCylwOB+FtkJsRTaKbIU4bEUUFdmgyI6ldAJylH52JRxHqeTkcD5b7ZT8QhVwnkJAECn92vmuIMVnJpxJWYrnut6XS+eLXHgGsNkV5wrtBAdYsJpNxWst+XiakJJtXlifYAKKBGw4KBRFkTgoxIFNXMeN2CnEOa9QHNhdX97FX+AOEyVf6I4yryuaLyKYzRasZisWiwWryYrZYsXPbMVidj77ma34mayuz7e15GExWbCYnMuUnXY+O6cDXNNWkwWr2c/1vh8Ws7XMtLVk2mr2w2L2c85zPYvJfNljHOBv3+xg+ur9fPd4Xzo2D3bL4VEVSinnseF6dijnMVJ6nkJhd7ielR2Hw8FfZrRyeyzuSBTl/aXLVvSqUsZZUGQKMAUgICqAe2XTpUuUai9mMVnwN/tf9PAz+xFgDij5Mg0uPd8ScMn7/mXeL28dZctYTVZM4rv3AoxLiORfK/ayMDmDxwZ3cP8Gis+jFr80mTFjxez+LV01u8NekmguehSVTjyu56L8i8oUOYpKDlyFcj4r57ODUtOl3i+ethdPK4UDx0XLXlTW9X7KwVMU2u0ktG5Y4brL3TYOLGIp9ePlwo+ZemZ/Gpislf8oKfPjp/iHSXnL+Jn9nMlBLLXuetPUfu2Yu/4Q/1i6m3fujK+x7YoIFnF+RVux1th2y+OORJEBlG4mGgGU7aa0KmUAUEpNB6YDdI7trOYMn1NhEvBz/ZrQrlxk4/r06dCUBcnpPDywfZ28L99sMlPfVJ/61vpGh1Khgyfz6Pf9Sp4c1omHBrQ3Opw6qVGgH/f1acOby/awOT2b2MiGRodU49zxk3gj0EFE2oiIHzAB+KJMmS+Au8SpF3DmctcnAIKsQcQ3j6db0250aNSBVg1a0SKwBQ0DGlLfWl8niWoanxjJ4ezzrNHjBXutRSkZmARu66HbThjp3uvb0DjQj9eWeuACgA+odqJQShUBDwPfATuABUqpbSIyVUSmuop9A+wH9gIzgAeru12t+oZ0aU7jQD/m6ZbaXsnuUCxKyaBvx1BahAQYHU6dFhxg5cH+7Vi9J4u1dfCHlVtOsiulvlFKdVRKtVNKveSa945S6h3XtFJKPeR6P1opleyO7WrV428xM7p7ON9vP0ZWboHR4Whl/LQ3i8wz+YyN1x0AeoPf9GpNy5AAXvlu14VbnOsI370aq7nFhKRIihyKT1P0YC3eZkFyOg3rWxncpZnRoWhAgNXMY4M7sDk9m6XbjxkdTo3SiaKOa98smITWjZi/Mb3O/UryZmfO2Vi6/Ri3xoXjb/Gme8Xqttt6RNA2NJDXvtuF3QONYb2VThQaE5JasT8rjw0HThkdiubyxebDFBY5dAeAXsZiNvH7IZ3YczyXJb8cNjqcGqMThcaN0S0I9rfojgK9yILkDK4Ja0C38BCjQ9HKuKFbC7qFN+CNZbspLHIYHU6N0IlCo76fhZHdW/L1lkzOnNOjehltR+ZZthw+w7gEXZvwRiaT8OSwzmScPl9nxnbRiUIDYEJiKwqKHHy+ue5Up73VwuQMrGZhZFy40aFoFejboSm92jbmXyv2kldQZHQ4HqcThQZAt/AQuoU34JMN+qK2kQqLHCxJPVzSxkXzTiLCH4d3Jiu3gPd+OmB0OB6nE4VWYnxiq5LTHpoxVuw8zqm8Qt12wgf0aNWIwdc0591V+8k+V7uHTNWJQisxMq4lAVYTn2zQF7WNsjA5nWbB/vTp0NToULQqeHJYJ3ILivhvLR8yVScKrUSDACs3Rbfki9TDdeK8q7c5npPPyt0nuC0+AotZH5q+oFOLYG6NC2fOT7V7yFT9adQuMjEpkrxCO1//etk+GzU3+2zTYewOxVjddsKnPDG4I3aH4q3le4wOxWN0otAuEt+6Ee2bBTFvY9247c9bKKVYkJxOfOtGtA0NMjoc7Qq0alKfiUmtmL8xnYMn84wOxyN0otAuIiJMSIxk06Fsdh/LMTqcOuOX9Gz2ncjTbSd81CMD22M1m/jH97uNDsUjdKLQLjGqezhWszBPX9SuMQuT06lnNXNTTEujQ9GuQrMGAdx9XRSfpx5h+5GzRofjdjpRaJdoEuTP0K4tWPxLBvk2++UX0KrlfKGdLzdncmN0GEH+ejAuX/VA33Y0CLDUysGNdKLQyjUhMZJsVw+mmmd9uy2T3IIixurTTj4tpL6Vqf3bsWLncZLTalcHmzpRaOW6rl1TIhrV06Pf1YAFGzNo1bg+Pds0NjoUrZru7t2G0GB//v5t7RrcqFqJQkQai8j3IrLH9dyonDKRIvKDiOwQkW0i8lh1tqnVDJNJGJ8Qydp9J2vtnRzeIP3UOdbtP8nY+AhExOhwtGqq52fm0YHt2ZB2ipW7TxgdjttUt0YxDViulOoALHe9LqsI+L1S6hqgF/CQiHSp5na1GjA2IRKToLsf96BFKRmIwG267UStMT6xFZGN6/Hqt7tw1JLBjaqbKEYC77um3wduLVtAKZWplNrkms4BdgC6W0wf0CIkgAGdmrEwJYMie93od78mORyKRSkZXN++KS0b1jM6HM1N/CzOwY22Z57l6y21o+FqdRNFc6VUJjgTAlDp4L4iEgV0B9ZXUmaKiCSLSPKJE7Wn6uarJiS14kROASt2Hjc6lFpn3f6THM4+z9gE3QFgbXNLbEs6twjm9aW7sNWCH1mXTRQiskxEtpbzGHklGxKRIOBT4HGlVIU3GiulpiulEpRSCaGhoVeyCc0DBnQKpVmwvz795AELktNpEGBhaJfmRoeiuZnJJPxhaCfSTp5jYXKG0eFU22UThVJqsFKqWzmPz4FjIhIG4Hou92eniFhxJom5SqnF7twBzbMsZhNjEyL4YddxMs+cNzqcWuPMeRvfbj3KyLhwAqxmo8PRPGDQNc3o0aoh/1y+2+fbI1X31NMXwCTX9CTg87IFxHkrxyxgh1LqH9XcnmaAcQmROBQsqgW/jLzFV78eoaDIodtO1GLFgxsdO1vAB+vSjA6nWqqbKF4GhojIHmCI6zUi0lJEvnGVuQ64ExgoIqmux43V3K5Wg1o3CeS69k2Yn5xea+7iMNqC5Aw6NQ8mOjzE6FA0D+rVtgl9O4byn5X7OJvvu+PRVytRKKVOKqUGKaU6uJ5PueYfUUrd6Jpeo5QSpVSMUirO9fim8jVr3mZ8YisyTp/np31ZRofi83Yfy2FzejZjE3Tbibrgj8M6kX3OxsxV+40O5arpltlalQzr2pyG9a26o0A3WJicjsUkjOqu7xKvC7qFh3BTdBgz1xwgK7fA6HCuik4UWpX4W8yM7h7B0u1HOemjH3ZvYLM7+OyXwwy6phlNgvyNDkerIb8b2pGCIgf/XrHX6FCuik4UWpVNSIrEZlcs3nTY6FB81spdJ8jKLWRsvG47UZe0Cw1ibHwEH68/RMbpc0aHc8V0otCqrGPzYHq0asi8jYdqVYdnNWlBcjpNg/zp30m3EaprHhvcAQTeXOZ7Q6bqzu+1KzIhqRV/XPQryQdPkxilezu9Elm5Bfyw8zj3Xt8Gi9lzv9FsNhsZGRnk5+d7bBva1flgdDi5+UX8unUb1mp+BgICAoiIiMBqtbopuorpRKFdkRExYbzw5XbmbUjXieIKLfnlMEUO5fG2ExkZGQQHBxMVFaXvqvIyRXYHu47mEBRgoXWTwKtej1KKkydPkpGRQZs2bdwYYfn0qSftitT3s3BLXEu+3nKEM+d9977wmqaUYkFyOnGRDWnfLNij28rPz6dJkyY6SXghi9lE02B/zpy3ca6w6KrXIyI0adKkxmqNOlFoV2xCYiT5NgdfbD5idCg+49eMM+w+lsu4GuoAUCcJ79U0yB+LycTRM9X7kq/J/7FOFNoViw4PoUtYAz363RVYkJxOgNXEiNgwo0PRDGY2CaHB/uQWFJHrI621daLQrpiIMCEpkm1HzrL18Bmjw/F6+TY7X2w+wg3dwmgQ4PkLj97gpZdeomvXrsTExBAXF8f69RWOLFChiRMnEhMTwxtvvOGBCCs2efJkFi1a5NZ1zpkzhyNHLtTApz3+IAf37ebo2QKfuINQX8zWrsrIuHBe+noHn2w4xEujoo0Ox6t9t+0oOflFjK0jo9itW7eOr776ik2bNuHv709WVhaFhYVVXr6oqIisrCzWrl3LwYMHPRhpzZkzZw7dunWjZcuWAMyaNYtTeYVknD7H2fwiQup59w8IXaPQrkpIPSs3RYfxReqRal2UqwsWJmcQ0agevdo2MTqUGpGZmUnTpk3x93e2PG/atGnJF2RUVBRZWc7+wpKTk+nfvz8Azz//PFOmTGHo0KHcddddDB06lOPHjxMXF8fq1auZMWMGiYmJxMbGctttt3HunLPR2rFjxxg1ahSxsbHExsaydu1aAD766COSkpKIi4vjgQcewG6/tJvvlJQU+vXrR3x8PMOGDSMz89LR6Mors2PHDpKSkkrKpKWlERMTA8ALL7xAYmIi3bp1Y8qUKSilWLRoEcnJydxxxx3ExcVx/vx5+vfvz77tm/G3mJk150Oio6Pp1q0bTz31VMl6g4KCeOaZZ4iNjaVXr14cO3asuv+aq6ZrFNpVm5DUisW/HObrXzP1KG0VyDh9jp/2ZfHYoA6YTDV/gfnPX25j+5EKxwm7Kl1aNuC5m7tW+P7QoUN54YUX6NixI4MHD2b8+PH069fvsutNSUlhzZo11KtXj7S0NEaMGEFqaqpzm126cP/99wPwpz/9iVmzZvHII4/w6KOP0q9fPz777DPsdju5ubns2LGD+fPn89NPP2G1WnnwwQeZO3cud911V8m2bDYbjzzyCJ9//jmhoaHMnz+fZ555htmzZ1epTGFhIfv376dt27bMnz+fcePGAfDwww/z7LPPAnDnnXfy1VdfMWbMGP7973/z2muvkZCQULJ+EUHlneTVF5/lx5/W0zaiOUOHDmXJkiXceuut5OXl0atXL1566SX++Mc/MmPGDP70pz9V/R/lRrpGoV21xKhGtA0N1KPfVeLTlMMoBbf1qBunncD5SzglJYXp06cTGhrK+PHjmTNnzmWXu+WWW6hXr/yxw7du3UqfPn2Ijo5m7ty5bNu2DYAVK1bw29/+FgCz2UxISAjLly8nJSWFxMRE4uLiWL58Ofv3X9xz665du9i6dStDhgwhLi6OF198kYyMjCqXGTduHAsWLABg/vz5jB8/HoAffviBnj17Eh0dzYoVK0rirMjOLan07N0Hh38wJrOZO+64g1WrVgHg5+fHiBEjAIiPjyctLe2yf0NP0TUK7aqJCBMSI/nrNzvZcyyHDs092z7A1zgcikWb0rmufRMiG9c3JIbKfvl7ktlspn///vTv35/o6Gjef/99Jk+ejMViweFwjiFdtg1AYGDFDdAmT57MkiVLiI2NZc6cOaxcubLCskopJk2axN/+9rdKy3Tt2pV169ZdVZnx48czduxYRo8ejYjQoUMH8vPzefDBB0lOTiYyMpLnn3++Su0c6vmZKbQ7OJV38XUcq9Vacgus2WymqMi4U7y6RqFVy+geEVjNomsV5Vh/4BTpp87XuQ4Ad+3axZ49F/ozSk1NpXXr1oDzGkVKSgoAn376aZXXmZOTQ1hYGDabjblz55bMHzRoEP/9738BsNvtnD17lkGDBrFo0SKOH3eOzHzq1KlLLop36tSJEydOlCQBm812ya//ysq0a9cOs9nMX/7yl5LaRHFSaNq0Kbm5uRfdORUcHExOTs4l+9WzZ0/WrllNYW42mafP8fEnn1TpNF1N04lCq5amQf4M6dKcTzdlUFDk2+MCu9vC5HSC/S0M69rC6FBqVG5uLpMmTaJLly7ExMSwfft2nn/+eQCee+45HnvsMfr06YPZXPWxwv/yl7/Qs2dPhgwZQufOnUvm//Of/+SHH34gOjqa+Ph4tm3bRpcuXXjxxRcZOnQoMTExDBky5JIL1X5+fixatIinnnqK2NhY4uLiSi6EV7XM+PHj+eijj0quTzRs2JD777+f6Ohobr31VhITE0vKTp48malTp5ZczC4WFhbG3/72NyaNGcGowb3p3DWGkSNHVvnvUlOkOvfwikhjYD4QBaQB45RSpysoawaSgcNKqRFVWX9CQoJKTk6+6vi0mvHj7hNMmr2Bf9/enRExLY0Oxyvk5NtIfGkZo3tE8Ncavn14x44dXHPNNTW6Ta360rLyyCsoolOL4Cp3Glne/1pEUpRSCRUsclWqW6OYBixXSnUAlrteV+QxYEc1t6d5oT7tmxLesJ4e/a6Ur3/NJN/mqDNtJ7TqaxESgF0pTnjhwGDVTRQjgfdd0+8Dt5ZXSEQigJuAmdXcnuaFTCZhXEIka/ZmkX7K9wZl8YQFyem0bxZEXGRDo0PRfESA1Uyj+n6czC3EVuQwOpyLVDdRNFdKZQK4nptVUO5N4I/AZfdeRKaISLKIJJ84caKa4Wk1ZWxCBCZBX9QG9h7PZdOhbMYlROjO+bQr0ryBPwo4luNdY4lcNlGIyDIR2VrOo0pXXERkBHBcKZVSlfJKqelKqQSlVEJoqB4FzFe0bFiPfh1DWZiSTpHdu34N1bSFKemYTcKt3cONDkXzMX4WM40D/TidZ6PA5j03h1w2USilBiulupXz+Bw4JiJhAK7n4+Ws4jrgFhFJA+YBA0XkIzfug+YlJiS14tjZAlbuqrs1wSK7g8WbDjOgUzOaBQcYHY7mg5oF+yMCx856z7WK6p56+gKY5JqeBHxetoBS6mmlVIRSKgqYAKxQSv2mmtvVvNDAzs1oGuTPvDp8+mnVnhOcyCnw+Ch2Wu1lNZtoGuRH9vlCzntJP2rVTRQvA0NEZA8wxPUaEWkpIt9UNzjNt1jNJsYmRPDDruMcO+td51hryoKNGTQJ9GNg54ou19V+aWlpdOvW7aJ5zz//PK+99lq113216zly5AhjxowBnA0Av/nGu7+emgb7YzYJR72kVlGtRKGUOqmUGqSU6uB6PuWaf0QpdWM55VdWtQ2F5pvGJ0RidygWpWRcvnAtczK3gOU7jzGqezjWKt4Hr9WMli1blrSU9oVEYTGZCA32JyffRl6B8bUK/WnW3CqqaSDXtm3CvI2HcDi8f0AWd1qSegSbXemedKvozJkzREVFlfT9dO7cOSIjI7HZbOzbt4/hw4cTHx9Pnz592Llz5yXLp6am0qtXL2JiYhg1ahSnTzvb+u7du5fBgwcTGxtLjx492LdvX0ktp7CwkGeffZb58+cTFxfH/Pnz6dChA8V3WDocDtq3b1/SFbqRmgb6YzU7h0w1enAj3Smg5nYTkiJ5bF4q6/af5Lr2TY0Op0YopViYnE5MRAidWnhR54j/mwZHt7h3nS2i4YaXq72akJAQYmNj+fHHHxkwYABffvklw4YNw2q1MmXKFN555x06dOjA+vXrefDBB1mxYsVFy991113861//ol+/fjz77LP8+c9/5s033+SOO+5g2rRpjBo1ivz8fBwOR0m/T35+frzwwgskJyfz73//G4CdO3cyd+5cHn/8cZYtW0ZsbCxNmxr/uTWZhGbB/hzOPk9OfhENDBzcSNcoNLcb1rUFIfWsfFKHxtTeduQsO4/m6NoEVNh2pLz548ePZ/78+QDMmzeP8ePHk5uby9q1axk7dmzJwENl+2o6c+YM2dnZJR3oTZo0iVWrVpGTk8Phw4cZNWoUAAEBAdSvX3nPvffccw8ffPABALNnz+buu+++sh32oEaBfvhZTBw9a2ytQtcoNLcLsJoZ1T2cj9cf4lReIY0D/YwOyeMWJKfjZzFxi7f1deWGX/5XqkmTJiWngYqdOnWKNm3aXFL2lltu4emnn+bUqVOkpKQwcOBA8vLyaNiwYcmgRVfiar5MIyMjad68OStWrGD9+vUX9U5rNJMILRoEcOjUOc6ct9GwvjHHkq5RaB4xISmSQruDxZtq/0XtfJudz1OPMLxrC0Lqe/fYxzUhKCiIsLAwli9fDjiTxLfffsv1119fbtmkpCQee+wxRowYgdlspkGDBrRp04aFCxcCzi//zZs3X7RcSEgIjRo1YvXq1QB8+OGH9OvXjwYNGhAREcGSJUsAKCgoKBk2tVh5XX7fd999/OY3v2HcuHFX1KttTQipZyXAaubo2XwcBtUqdKLQPKJziwbERTZk/sZ0wy/EedqyHcc4c96m206U8sEHH/Diiy8SFxfHwIEDee6552jXrl25ZYu76y4e1wFg7ty5zJo1i9jYWLp27crnn1/SRIv333+fJ598kpiYGFJTU0uGIP3www956623iImJoXfv3hw9evSi5QYMGMD27dtLLmaDs2aTm5vrVaediomrVlFY5OB0mcGNaiwGbz6IdTfjvm3ehkNMW7yFT397LfGtGxsdjsfcNXsDe4/lsPqpgZgNGBe7LN3N+JVLTk7miSeeKKmheBulFPtO5GGzO+jUPLhk/HVf6WZc0yp0c2xLAv3Mtbr78SPZ51m95wRj4iO8IkloV+7ll1/mtttuq3ToVKOJCC1CArDZHZzMq/lGeDpRaB4T6G/h5tiWfPVrJjn5NqPD8YjFmzJQCsbUseFOa5Np06Zx8ODBcq+heJMgfwvBAVaO5xRQ5KjZjjd1otA8akJSK87b7Hyx+YjRobidUs4W6L3aNqZVk8pvwdQ0d2jRwB+7Q5GVU7PXKnSi0DwqNiKEzi2Ca+Xpp41pp0k7eY6xujah1ZB6fhYa1rOSlVuArQa789eJQvMoEWFCYiRbDp9h6+EzRofjVguS0wnyt3BDdAujQ9HqkOYNAlAKTuTU3LUKnSg0j7u1ezh+FlOtGv0ut6CIb7ZkMiImjPp+ut2qVnP8rWYaBVo5mVdYY9cqdKLQPK5hfT9u7NaCJamHOV/oPaN2Vcc3v2ZyrtCu206Uo3///nz33XcXzXvzzTd58MEHq73uoKCgq1ru2WefZdmyZSWxlG2E52uKB8XKOV8zPcvqRKHViPGJrcjJd/4Krw0WpqTTNjSQHq0aGR2K15k4cSLz5s27aN68efOYOHGiQRHBCy+8wODBg4HakSj8LCaaBvpxrtDO7mM5l1+gmnSi0GpEr7aNiWpSv1acftp/IpeNaacZGx9ZYQd4ddmYMWP46quvKChwnkNPS0vjyJEjl9x++tRTT/Gf//yn5PXzzz/P66+/DsCrr75KYmIiMTExPPfcc5dsQynFk08+Sbdu3YiOji5pYQ3w97//nejoaGJjY5k2bRoAkydPZtGiRbz11lscOXKEAQMGMGDAAGbNmsUTTzxRsuyMGTP43e9+574/hgeFuoZMfe27XR7fVrVOropIY2A+EAWkAeOUUqfLKdcQmAl0AxRwj1JqXXW2rfkWEWF8Yite+XYne4/n0r7Z1Z1CMJrDofhg3UFMAqN7hBsdzmW9suEVdp66dCyH6ujcuDNPJT1V4ftNmjQhKSmJb7/9lpEjR5b0Cls2qU6YMIHHH3+85JTUggUL+Pbbb1m6dCl79uxhw4YNKKW45ZZbWLVqFX379i1ZdvHixaSmprJ582aysrJITEykb9++pKamsmTJEtavX0/9+vU5derURdt89NFH+cc//sEPP/xA06ZNycvLIyYmhr///e9YrVbee+893n33XTf+tTzHYjYRFGBh6fZ0fjl0mu4erN1Wt0YxDViulOoALHe9Ls8/gW+VUp2BWGBHNber+aDb4sOxmIQFyb5Xq8i32Zm34RBD3viROWvTuDE6jOYNAowOy2uVPv1U0Wmn7t27c/z4cY4cOcLmzZtp1KgRrVq1YunSpSxdupTu3bvTo0cPdu7cyZ49ey5ads2aNUycOBGz2Uzz5s3p168fGzduZNmyZdx9990lXYs3blx51zGBgYEMHDiQr776ip07d2Kz2YiOjnbTX8HzgvwtNAn041UP1yqqe7vGSKC/a/p9YCVw0U8NEWkA9AUmAyilCgFjerbSDNUsOIBB1zTj05QM/jC0E34W7z/zeTqvkLnrDzJn7UGycgvoEtaAf06I48boMKNDq5LKfvl70q233srvfvc7Nm3axPnz5+nRo0e55caMGcOiRYs4evQoEyZMAJynlZ5++mkeeOCBCtdfUR91SqkrPh1433338de//pXOnTt7ZaeAlTGJ8NCA9rzw1XbW7Mni+g6eGXCpukdqc6VUJoDrubwR5dsCJ4D3ROQXEZkpIoHV3K7moyYkteJkXiHLdhwzOpRKHTp5juc+30rvl1fw2tLddG3ZgI/v68nXj17PyDg9JvblBAUF0b9/f+65555KL2JPmDCBefPmsWjRIsaMGQPAsGHDmD17Nrm5uQAcPny4ZIS6Yn379mX+/PnY7XZOnDjBqlWrSEpKYujQocyePbvkYnXZU09waTfjPXv2JD09nY8//tjQC+5X645erQhvWI9Xv9vpsZ6aL1ujEJFlQHktip65gm30AB5RSq0XkX/iPEX1/yrY3hRgCkCrVq2quAnNV/TtEErLkAA+2XDIK3+Vp6ZnM2PVfv63NROzSRgZF879fdp61/CmPmLixImMHj36kjugSuvatSs5OTmEh4cTFub8PAwdOpQdO3Zw7bXXAs6k89FHH9Gs2YXfoaNGjWLdunXExsYiIvz973+nRYsWDB8+nNTUVBISEvDz8+PGG2/kr3/960XbnDJlCjfccANhYWH88MMPAIwbN47U1FQaNfK9u9j8LWYeG9yBPy76le+2Hb38AlehWt2Mi8guoL9SKlNEwoCVSqlOZcq0AH5WSkW5XvcBpimlbrrc+nU347XTG9/v5q0Ve1j15AAiGxvfR5LDoVi+8zgzVu1nQ9opggMs/KZXayb3jvLJ6xC6m/ErN2LECJ544gkGDRpkdChXpPh/XWR3MOzNVYgIy3/f3+u6Gf8CmOSangRcMrqIUuookC4ixQlkELC9mtvVfFhxI7WFBl/UzrfZ+WTDIQa/8SP3f5DM4ezz/L8RXVj39CCeGt7ZJ5OEdmWys7Pp2LEj9erV87kkUZrFbOIPQzux93iuZ9ZfzeVfBhaIyL3AIWAsgIi0BGYqpW50lXsEmCsifsB+wLeuGGluFdGoPn07hLIgOYPHBnes8XEcTucV8uHPB/lgXRpZuYV0C2/AWxO7c2O3Flj0tYc6pWHDhuzevdvoMNxieLcWxESEcNAD665WolBKncRZQyg7/whwY6nXqYBbq0Kab5uQGMlv527ix93HGdi5eY1s8+DJPGatOcCC5HTybQ4GdArl/r5tubZtE91wTvN5IsIfh3XmSw+sW/dmphli0DXNaRrkx7wN6R5PFL8cOs30Vfv5dttRLCbh1rhw7u/blo7N9QVqrXbx1O2xOlFohvCzmLitRwQz1xzg+Nl8mrn5eoDDoVi24xgzVu9nY9ppGgRY+G2/dkzuHeX2bWlabacThWaY8YmRvLtqP4s2ZfBg//ZuWWe+zc7iTYeZuXo/+7PyCG9Yj2dHdGF8YiSB/vrjrmlXQ1+50wzTNjSIpDaNmb8xHYejeg2FTuUV8s9le7ju5RX832dbCPS38K+J3fnxyf7cc30bnSRqWHFL63bt2tGlSxduvPHGGr9ovHLlStauXVuj2yxPWloa3bp1MzqMatFHj2aoiUmRPDF/Mz8fOEnvdld+fjUty3mBemGK8wL1wM7NuL9PW3q1bawvUBtEKcWoUaOYNGlSSWO71NRUjh07RseOHau0DrvdjtlsrvB1VaxcuZKgoCB69+59Rctpl9I1Cs1QN3QLo0GA5YrH1N506DRTP0xhwOsrmb8xnZGx4Xz/RF9mT07k2nb6LiYj/fDDD1itVqZOnVoyLy4ujj59+rBy5UpGjBhRMv/hhx9mzpw5AERFRfHCCy9w/fXXs3DhwkteL126lGuvvZYePXowduzYki4+oqKieO655+jRowfR0dHs3LmTtLQ03nnnHd544w3i4uJYvXr1RTH++OOPxMXFERcXR/fu3cnJyWHlypX07duXUaNG0aVLF6ZOnYrDNYJcRdtOSUmhX79+xMfHM2zYMDIzM0vmx8bGcu211/L222977G9dU3SNQjNUgNXMqO7hfLIhndN5hTQK9KuwrMOh+H7HMWas2k/ywdOE1LPyYP92TOodVTLil3axo3/9KwU73NvNuP81nWnxf/9X4ftbt24lPj7+qtYdEBDAmjVrAJg2bVrJ66ysLEaPHs2yZcsIDAzklVde4R//+AfPPvssAE2bNmXTpk385z//4bXXXmPmzJlMnTqVoKAg/vCHP1yynddee423336b6667jtzcXAICnJ+fDRs2sH37dlq3bs3w4cNZvHgx/fv358UXX7xk208//TSPPPIIn3/+OaGhocyfP59nnnmG2bNnc/fdd/Ovf/2Lfv368eSTT17V38Kb6EShGW58YiveX3eQJamHufu6Npe8n2+z8+mmDGauPsCBrDwiGtXj+Zu7MDZBX6CubcaPH1/u659//pnt27dz3XXXAVBYWFjSFxTA6NGjAYiPj2fx4sWX3c51113H7373O+644w5Gjx5NRISzt4CkpCTatm0LOPuqWrNmDQEBAeVue9euXWzdupUhQ4YAztNjYWFhnDlzhuzsbPr16wfAnXfeyf/+97+r/pt4A32UaYbr0rIBsREhzNuQzuTeUSWnjU7mFrhaUB/kVF4hMREh/Pv27gzvqltQV1Vlv/w9pWvXrixatKjc9ywWS8npHID8/PyL3g8MDCz3tVKKIUOG8Mknn5S7Xn9/fwDMZjNFRZcfR3ratGncdNNNfPPNN/Tq1atkPO2ypyxFpMJtb9myha5du7Ju3cVjsGVnZ9e6U5/6aNO8wvjEVuw6lkNqejYHsvL405It9H55BW8u20OPVg2ZP6UXnz90HSNiWuok4eUGDhxIQUEBM2bMKJm3ceNGfvzxR1q3bs327dspKCjgzJkzLF++vErr7NWrFz/99BN79+4F4Ny5c5e9i6psd+Kl7du3j+joaJ566ikSEhLYudN5em7Dhg0cOHAAh8PB/Pnzuf766yvcdqdOnThx4kRJorDZbGzbto2GDRsSEhJScgpt7ty5VdpHb6aPOM0r3BLXkvp+Zn770SYGvr6SBRszGNU9nGW/68vMSYn01N1s+AwR4bPPPuP777+nXbt2dO3aleeff56WLVsSGRnJuHHjiImJ4Y477qB79+5VWmdoaChz5sxh4sSJxMTE0KtXr5Iv94rcfPPNfPbZZ+VezH7zzTfp1q0bsbGx1KtXjxtuuAGAa6+9lmnTptGtWzfatGnDqFGjKty2n58fixYt4qmnniI2Npa4uLiS23Hfe+89HnroIa699lrq1at3FX9F71KtbsY9TXczXrf8+cttLN50mDt7teau3q31BeqrpLsZvzorV67ktdde46uvvjI6lCor738tIm7vZlxfo9C8xrMjuvDsiC665qBpXkYnCs1r6AShGal///7079/f6DC8kr5GoWm1kDefUtbcoyb/xzpRaFotExAQwMmTJ3WyqMWUUpw8ebKkoaCn6VNPmlbLREREkJGRwYkTJ4wORfOggICAkoaCnlatRCEijYH5QBSQBoxTSp0up9wTwH2AArYAdyul8suW0zSt+qxWK23aXNrCXdOuVnVPPU0DliulOgDLXa8vIiLhwKNAglKqG2AGJlRzu5qmaVoNqW6iGAm875p+H7i1gnIWoJ6IWID6wJFqblfTNE2rIdVNFM2VUpkArudmZQsopQ4DrwGHgEzgjFJqaTW3q2maptWQy16jEJFlQIty3nqmKhsQkUY4ax5tgGxgoYj8Rin1UQXlpwBTXC8LRGRrVbbjQSHAGYNjcAe9H96ntuxLbdmPpkCW0UG4QSd3r/CyiUIpNbii90TkmIiEKaUyRSQMOF5OscHAAaXUCdcyi4HeQLmJQik1HZjuKpvs7qboV0pEpiulply+pHfT++F9asu+1KL9MPz7xh1ExO39HlX31NMXwCTX9CTg83LKHAJ6iUh9cTa9HQTsqOZ2a9KXRgfgJno/vE9t2Zfash9aBarVKaCINAEWAK1wJoSxSqlTItISmKmUutFV7s/AeKAI+AW4TylVUIX114oMr2ma96st3zee2A+v7j1WRKa4TkVpmqZ5VG35vvHEfnh1otA0TdOMV6f6ehKR4SKyS0T2isg017xYEVknIltE5EsRaVDVZV3zG4vI9yKyx/XcqAb2Y7aIHC99R5iIvCoiO0XkVxH5TEQa+uh+PC8ih0Uk1fW40Uf3I05EfnbtQ7KIJPnAfkSKyA8iskNEtonIY675Y12vHSJS4SkNb9oXzc2UUnXigbNF+D6gLeAHbAa6ABuBfq4y9wB/qeqyrvf+DkxzTU8DXqmBfekL9AC2lpo3FLC4pl8pLw4f2Y/ngT9czf/Sy/ZjKXCDa/pGYKUP7EcY0MM1HQzsdh0j1+C85XIlzh4WvP5/oh/ufdSlGkUSsFcptV8pVQjMw9m+oxOwylXme+C2K1gWqt463W2UUquAU2XmLVVKFY8q/zNQXm9hXr8fVeQL+6GA4tppCOX3RuBt+5GplNrkms7BeXdiuFJqh1Jq12UW96p9qeDsQZVqNt5WM6pgX/7iOnuQKiJLXTcQeWxf6lKiCAfSS73OcM3bCtzimjcWiAQQkZYi8s1lloUqtE43wD3A/8Bn9+Nh10Ewu/gD7IP78Tjwqoik4+yZ4Gnwnf0QkSigO7C+kjJeuS8iYgbeBm7AWSOaKCJdqFrfdBUtS1WWr8F9eVUpFaOUigO+Ap715L7UpURR3vBpCueX6kMikoKzul0IoJQ6oly391ayrNcRkWdw3oY8F3xyP/4LtAPicHb58jr45H78FnhCKRUJPAHMAt/YDxEJAj4FHldKna2onBfvS0W1m6rUbLyqZlRRPGX+L4GU/7d2277UpUSRgau24BIBHFFK7VRKDVVKxQOf4DzPWqVlXdPHxNkqHam4dXqNEJFJwAjgDqVUeR8cr98PpdQxpZRdKeUAZuD8sJfl9fuBswHqYtf0QnxkP0TEijNJzFVKLb5c+VK8aV8qqt2UW7Px1prR5eIRkZdcNdY7cNUoPLUvdSlRbAQ6iEgbEfHD2dX5FyJS/GExAX8C3qnqsq73qtI63eNEZDjwFHCLUupcBcV8YT/CSr0chfPUYFlevx84vyT7uaYHAnvKKeNV+yEigrPms0Mp9Y8rXNyb9uWKajdeXDOCSuJRSj3jqrHOBR52zfPMvhh9Nb0mHzjvPtmNs9bwjGveY655u4GXudC2pCXwTWXLuuY3wXmOb4/ruXEN7McnOE/L2HD+SrgX2Ivz10Oq6/GOj+7HhzgHt/oV5xdMmI/ux/VACs67f9YD8T6wH9fj/CL5tdTn6EacCTsDKACOAd95874A1xbH6Hr9tOuxq9TnKQzYVdVlXdOXXb6m9qVMmdaUuuPOE/uiG9xpmlariHPcm904+5U7jLO2czswGTiplHrZdQdQY6XUH6uyrFJqm4i8ernla3BfCpVSe1xlHsF5i/8YT+2LThSaptU64myo+SbO9h2zlVIvSdX7prtkWdf8cpc3aF8+xXlrvwM4CExVSh321L7oRKFpmqZVqi5dzNY0TdOugk4UmqZpWqV0otA0rVaooKuL+XKhg8k0EUmtZPnuIqJEZFiNBe0jLjsUqqZpmrcr1V3FEJy38m4UkS+UUuNLlXmdysf2ngiscT1/58FwfY6uUWiaVhtU1l1FcWPCcTjbvFzC9f4YnLfQDhWRANf8KLm4+/g/iMjzrulEV59k68TZzX95jUNrBZ0oNE2rDSrrrgKgD3CsuO1BOa4DDiil9uHsTr3ccVDKeA/nbanXAvYrjtiH6EShaVptcLnuKiZSQW2i1PvzXNPzXK8r3phzYLBgpdRa16yPqxamb9LXKDRNqw0q7JTQ1UJ5NBBf/KaIvIezG/UjwM04x6G5xdX7sgBNRCQYZ0/MpX9QBxSvwjO74Z10jULTtNqgsk4JBwM7lVIZxYWVUncrpeJcLZgHA5uVUpFKqSilVGucPejeirNvq2Yi0kRE/HH2zoxS6jSQIyK9XKucUAP7aBidKDRN83nKObrjwzjvVtoBLFBKbXO9PYHLn3b6rMy8T3H2i2QDXsDZseNXwM5SZe4FpovIOpw1jMruqPJpugsPTdO0qyAiQUqpXNf0NJy9sT5mcFgeoa9RaJqmXZ2bRORpnN+jB3HeWlsr6RqFpmmaVil9jULTNE2rlE4UmqZpWqV0otA0TdMqpROFpmmaVimdKDRN07RK6UShaZqmVer/A0VYfEo681ElAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "df.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other ways to subset data "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assume that we interested in these 3 points only"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "pt1 = (606200, 6905480)\n",
    "pt2 = (606300, 6905410)\n",
    "pt3 = (606400, 6905520)\n",
    "pts_x = [pt1[0], pt2[0], pt3[0]]\n",
    "pts_y = [pt1[1], pt2[1], pt3[1]]\n",
    "elem_ids = dfs.find_nearest_elements(pts_x, pts_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use these element ids either when we select the data from the complete dataset using the method isel() as shown above or already when we read the data from file (particular useful for files larger than memory)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|█████████████████████████████████████████████████████████████████████████████████| 9/9 [00:00<00:00, 13878.21it/s]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<mikeio.Dataset>\n",
       "Dimensions: (9, 3)\n",
       "Time: 1985-08-06 07:00:00 - 1985-08-07 03:00:00\n",
       "Items:\n",
       "  0:  Surface elevation <Surface Elevation> (meter)\n",
       "  1:  U velocity <u velocity component> (meter per sec)\n",
       "  2:  V velocity <v velocity component> (meter per sec)\n",
       "  3:  Current speed <Current Speed> (meter per sec)"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds_pts = dfs.read(elements=elem_ids)\n",
    "ds_pts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Select area\n",
    "Let's take the area North of y=6905480"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "yc = dfs.element_coordinates[:,1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "elem_ids = dfs.element_ids[yc>6905480] "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And find the maximum average current speed in this area in the last time step"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "U velocity <u velocity component> (meter per sec)\n",
      "Max current speed in area is found in [ 606674.6 6905718.3] and is 0.039m/s\n"
     ]
    }
   ],
   "source": [
    "item_num = 1\n",
    "print(ds.items[item_num])\n",
    "subset = ds.data[1][:,elem_ids] \n",
    "subset_timeavg = subset.mean(axis=0) \n",
    "idx = subset_timeavg.argmax()\n",
    "coords = dfs.element_coordinates[idx,0:2].round(1)\n",
    "print(f'Max current speed in area is found in {coords} and is {subset_timeavg[idx]:.3f}m/s')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us save the time averaged subset to a dfsu file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\users\\jem\\source\\mikeio\\mikeio\\dfsu.py:2342: UserWarning: No start time supplied. Using start time from source: 1985-08-06 07:00:00 as start time.\n",
      "  warnings.warn(\n",
      "100%|██████████████████████████████████████████████████████████████████████████████████| 1/1 [00:00<00:00, 8355.19it/s]\n"
     ]
    }
   ],
   "source": [
    "outfilename1 = \"HD2D_north.dfsu\"\n",
    "data = []\n",
    "data.append(subset_timeavg.reshape(1,-1))\n",
    "items = ds.items[item_num]\n",
    "dfs.write(outfilename1, data, items=[items], elements=elem_ids)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Create a new dfsu file\n",
    "\n",
    "* Subset of items\n",
    "* Subset of timesteps\n",
    "* Renamed variables\n",
    "\n",
    "First inspect the source file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████████████████████████████████████████████████████████████████████████████| 9/9 [00:00<00:00, 7420.63it/s]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<mikeio.Dataset>\n",
       "Dimensions: (9, 884)\n",
       "Time: 1985-08-06 07:00:00 - 1985-08-07 03:00:00\n",
       "Items:\n",
       "  0:  Surface elevation <Surface Elevation> (meter)\n",
       "  1:  U velocity <u velocity component> (meter per sec)\n",
       "  2:  V velocity <v velocity component> (meter per sec)\n",
       "  3:  Current speed <Current Speed> (meter per sec)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds = dfs.read()\n",
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████████████████████████████████████████████████████████████████████████████| 5/5 [00:00<00:00, 9991.20it/s]\n"
     ]
    }
   ],
   "source": [
    "from mikeio.eum import ItemInfo, EUMType\n",
    "from mikeio import Dataset\n",
    "\n",
    "sourcefilename = filename\n",
    "outfilename2 = \"HD2D_selected.dfsu\"\n",
    "\n",
    "starttimestep = 4\n",
    "time = ds.time[starttimestep:]\n",
    "\n",
    "data = []\n",
    "data.append(ds['U velocity'][starttimestep:,:])\n",
    "data.append(ds['V velocity'][starttimestep:,:])\n",
    "\n",
    "items = [ItemInfo(\"eastward_sea_water_velocity\", EUMType.u_velocity_component),\n",
    "         ItemInfo(\"northward_sea_water_velocity\",EUMType.v_velocity_component)]\n",
    "\n",
    "newds = Dataset(data,time,items)\n",
    "\n",
    "dfs.write(outfilename2, newds) # Note, this method was previosly named create"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read the newly created file to verify the contents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████████████████████████████████████████████████████████████████████████████| 5/5 [00:00<00:00, 9972.19it/s]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<mikeio.Dataset>\n",
       "Dimensions: (5, 884)\n",
       "Time: 1985-08-06 17:00:00 - 1985-08-07 03:00:00\n",
       "Items:\n",
       "  0:  eastward_sea_water_velocity <u velocity component> (meter per sec)\n",
       "  1:  northward_sea_water_velocity <v velocity component> (meter per sec)"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "newdfs = Dfsu(outfilename2)\n",
    "\n",
    "newds2 = newdfs.read()\n",
    "newds2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Write mesh from dfsu file\n",
    "Don't you have the original mesh? No problem - you can re-create it from the dfsu file... "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "outmesh = 'mesh_from_HD2D.mesh'\n",
    "dfs.to_mesh(outmesh)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Clean up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "os.remove(outfilename1)\n",
    "os.remove(outfilename2)\n",
    "os.remove(outmesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "fa576ebcd40e010bdc0ae86b06ce09151f3424f9e9aed6893ff04f39a9299d89"
  },
  "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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}