{
 "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, 8998.51it/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/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAABPH0lEQVR4nO3dd3wUdf748dd7W0JIQoCEUBJI6FKSAEkAkSpNRThU2llARfRsp955euf91K+ed7bzvLPcCYLYjmpHPRAEAaUl9F4DKZRAICSBlN39/P7YTQiQhJLdzG7yeeo+dnb2MzPvCTv73s/MfD4fUUqhaZqmaZUxGR2Apmma5tt0otA0TdOqpBOFpmmaViWdKDRN07Qq6UShaZqmVUknCk3TNK1KFqMDqEx4eLiKiYkxOgxN0zS/kpqaelwpFeHJdfpsooiJiSElJcXoMDRN0/yKiBz09Dr1qSdN0zStSjpRaJqmaVXSiULTNE2rkk4UmqZpWpV0otA0TdOq5JFEISLDRWSXiOwVkacreL+liCwVkQ0isllEbvTEdjVN0zTvq3aiEBEz8A5wA9AJmCAinS4o9mdgrlKqGzAeeLe62/UVxXan0SFomqZ5lSdqFMnAXqXUfqVUMTAbGHVBGQWEuqcbAFke2K6hlFI8+9VWev1tCdl5RUaHo2ma5jWeSBQtgPRyrzPc88p7HrhDRDKA74BHPLBdQ726cBcfrTpITkExH/6SZnQ4mqZpXlNTF7MnADOVUlHAjcDHInLRtkVkioikiEhKdnZ2DYV25f7z0z7+vWwft/dsyfDOTfl49UEKiuxGh6VpmuYVnkgUmUB0uddR7nnl3QvMBVBKrQICgfALV6SUmqqUSlRKJUZEeLSrEo/575pDvPz9TkbGN+eFUV24v39rcs+WMGdd+qUX1jRN80OeSBTrgHYiEisiNlwXq7++oMwh4HoAEbkGV6Lw3SpDJb7elMUzX25hUMcm/H1sPGaT0K1lQ5JjGzF95QFKHPrCtqZptU+1E4VSyg48DCwEduC6u2mbiLwgIiPdxX4H3Ccim4BZwCSllKrutmvS0p3HeGLORpJiGvHu7d2xms/96e7v15rMU2f5dvNhAyPUNE3zDvHV7+vExETlK73Hrtl/grtmrKVdZDCz7utFSKD1vPedTsWwN5djMZv47tHrEBGDItU0ra4TkVSlVKIn16lbZl/C1sxcJn+YQlTDenx4d/JFSQLAZBKm9GvNjsOnWbHnuAFRapqmeY9OFFXYeyyfu2asJbSelU8m96RxcEClZUcltCAyNID3lu+rwQg1TdO8TyeKSmScPMOd09dgEuGTyT1p1qBeleVtFhP39Inl570n2JqZW0NRapqmeZ9OFBXIzivizulrKSiy89E9ycSG17+s5Sb0bElIgIX3lu/3coSapmk1RyeKC+SeLeGuGWs5klvIB3cn0al56KUXcgsNtPLrni35dnMW6TlnvBilpmlazdGJopwzxXbumbmOvcfymHpXD3q0anTF67i7TyxmkzB95QEvRKhpmlbzdKJwK7Y7eeCT9Ww4dJJ/je9G33ZX1zK8aYNARiW0YPa6Q+QUFHs4Sk3TtJqnEwXgcCoem7OB5buzefmWOG7o2qxa65vSrzWFJU4+XnXQQxFqmqYZp84nCqUUf/p8C99tOcKfb7qGsUnRl17oEtpHhnB9xyZ8uCqNwhKHB6LUNE0zTp1OFEop/vrdDuakpPPooLZM7tvaY+ue0q81OQXFzEvN8Ng6NU3TjFCnE8U7S/cybcUBJl0bw+ND2nt03cmxjUiIDuP9FftxOH2zmxRN07TLUWcTxUer0nh90W5u6daCZ0d08nj/TCLCA/1bc/DEGRZuO+LRdWuaptWkOpkovtiQwbNfbWNIp0hevS0Ok8k7nfgN6dSUmMZBvPfTPny180VN07RLqXOJ4oftR/n9vM1c26Yxb03ohsXsvT+B2STc1681mzJyWXMgx2vb0TRN86Y6lSh+2Xech/67ni4tGjD1rkQCrWavb/PW7lGEB9t47yfdWaCmaf6pziSKTemnuO/DFGIaBzFzUhLBAZYa2W6g1czE3jEs3ZXNriN5NbJNTdM0T6oTiWL30TwmfrCWxsEBfHxvTxrWt9Xo9u/s3Yp6VjNTdWeBmqb5IY8kChEZLiK7RGSviDxdSZmxIrJdRLaJyH89sd3LkZ7j6i7cZjbxyb09iQwNrKlNlwkLsjE+OZqvNmZyOPdsjW9f0zStOqqdKETEDLwD3AB0AiaISKcLyrQD/gj0UUp1Bh6r7nYvx7HThdz+/hpXdxr39qRl46Ca2GyF7r0uFgXM0J0FaprmZzxRo0gG9iql9iulioHZwKgLytwHvKOUOgmglDrmge1W6dSZYu6cvpYT+UV8eE8yHZqGeHuTVYpqGMSIuGbMWptO7tkSQ2PRNF+xat8JhrzxExsOnTQ6FK0KnkgULYD0cq8z3PPKaw+0F5GfRWS1iAyvaEUiMkVEUkQkJTs7+6oDKiiyM+mDdRw4UcC0uxJJiA676nV50pR+rckvsvPfNYeMDkXTDLdm/wnumbmOPcfy+et3O3RbIx9WUxezLUA7YAAwAZgmImEXFlJKTVVKJSqlEiMirq6b78ISB1M+TmFLZi5vT+jGtW3DqxG2Z3Vu3oC+7cKZ8fMBiuy6s0Ct7kpJy+HumetoHhbIY4PbsS7tJMt2X/2PQ827PJEoMoHyXa5GueeVlwF8rZQqUUodAHbjShweZXc4eXTWBn7ee4JXb41jaOemnt5Etd3frw3ZeUV8tSHL6FA0zRCpB08yccZamoYGMuu+Xjw4oC3Rjerx+sJdOHW/aD7JE4liHdBORGJFxAaMB76+oMyXuGoTiEg4rlNRHr1X1OlUPPXZFhZtP8rzN3fi1h5Rnly9x/Rp25jOzUN5b/k+fVBodc7G9FNMmrGWiJAA/ntfL5qEBmKzmHh8cHu2ZZ3mf7pfNJ9U7UShlLIDDwMLgR3AXKXUNhF5QURGuostBE6IyHZgKfCkUupEdbddLgZeWLCdz9Zn8MSQ9kzqE+upVXuciDClX2v2ZRewZKfXr+lrms/YkpHLndPX0LC+jVlTetG0wblb1UcltKBdk2D+vmgXdofTwCi1injkGoVS6julVHulVBul1Evuec8qpb52Tyul1BNKqU5Kqa5Kqdme2G6pNxfvYeYvadx7XSyPDGrryVV7xU1dm9EirB5Tl+tuPbS6YWtmLndMX0ODelZmTelFswb1znvfbBJ+N7QD+7IL+GLDhWeuNaP5fcvsGSsP8M8lexjTI4o/33SNx7sL9waL2cTkvrGsSztJ6kHdWaBWu23POs0d09cQHGBh1n29aBFWr8JywzpHEhfVgDcX79E3e/gYv04U81LSeWHBdm7o0pS/3dLVL5JEqXFJ0YQFWXnvJ92th1Z77TqSxx3T11DPambWfb2IblR5o1cR4clhHcg8dZbZa9MrLafVPL9NFP/bepinPttM33bhvDk+wavdhXtDkM3CXb1a8cOOo+zLzjc6HE3zuD1H8/j1tNVYzcKs+3pdVs8I17UNp1frRrz1417OFNtrIErtcvjXt6vbyj3HeXTWRhKiw3jvzh4EWLzfXbg33HVtDDazifdX6FqFVrvsPZbPhGlrMJlcSSImvP5lLVdaqzieX8TMX9K8G6R22fwuUaw/dJIpH6fQOqI+H0xKJshWM92Fe0N4cAC39Yjis9RMjuUVGh2OpnnE/ux8fj1tNQCz7utF64jgK1q+R6tGDOrYhPd+2q+7u/ERfpUodhw+zaQZa2kS4uouvEGQ1eiQqu2+vq0pcTr5UP960mqBtOMFTJi2GodT8d/7etK2yZUliVK/G9qe3LMlurbtI/wmUaQdL+DO6WsJsln4ZHJPIkICjA7JI2LC6zO8c1M+XnWQ/CJ9TlbzX4dOnGHCtNUU2518el9P2kdefUecnZs3YERcM6avPMDx/CIPRqldDb9IFEdyXd2FO5Xik8nJRDU0rrtwb5jSrzWnC+3MXqs7C9T8U3qOK0mcLXHw6eRedGwaWu11PjGkPUV2J+8u1e2NjObziSKnoJg7pq8h92wJH96dTNsmxnYX7g3dWjakZ2wjZqw8QIlular5mcxTZ/n1+6vJKyzhk3t70ql59ZMEQOuIYG7rHsUnqw+SeUoP+GUkn04UeYUlTJyxlvScM0yfmEjXqAZGh+Q19/dvTVZuIQs2684CNf9xOPcsE6au5tSZEj6Z3JMuLTx7jD462NV36FtL9nh0vdqV8dlEoRRM/jCFHYdP8+87utOzdWOjQ/KqAe2b0D4ymPd+2q/75df8wtHThUyYupqcgmI+uieZuKgwj2+jRVg9bu/VknmpGezX7Y0M47OJ4lBOAWvTcvj72HgGdYw0OhyvM5mEKf3asPNIHsv3HDc6HE2r0jF3ksjOK+LDe5Lo1rKh17b14IC22Mwm/rFY1yqM4rOJ4nShnRdHdWFUwoWD5dVeI+Ob0zQ0kPd+0hfvNN+VnVfEr99fw5HThcy8J5kerRp5dXsRIQHcc10M32zKYnvWaa9uS6uYzyaK1hH1uaNXK6PDqFE2i4l7rovhl30n2JKRa3Q4mnaRE/lF3P7+ajJPnuWDSUkkxXg3SZSa0rcNoYEW3vhhV41sTzufzyaK+n7c4ro6JiS3JCTAwnu6C3LNx+QUFHP7+2s4lHOG6ZMSa/S6YYMgK/f3b8PiHcdIPXiyxrarufhsoqirQgKt/LpXS77bcphDJ84YHY6mAXDqTDF3vL+GA8cLeP+uJK5tU/Nj0d/dJ4bwYBuvLdypb/ioYTpR+KB7+sRiNgnTV+ruCzTj5Z4p4Y7pa9h7LJ+pdyVyXbuaTxLg6nH54YFtWb0/h5/3emyATO0yeCRRiMhwEdklIntF5Okqyt0qIkpEEj2x3doqMjSQXyW0YE5KOjkFxUaHo9VhpwtLuGvGGnYdyeO9O3vQv32EofFM6NmSFmH1dK2ihlU7UYiIGXgHuAHoBEwQkU4VlAsBfgusqe4264Ip/VpTWOLko1VpRoei1VGlDV63Hz7Nv2/vwcCOTYwOiQCLmd8ObsemjFwWbT9qdDh1hidqFMnAXqXUfqVUMTAbGFVBuReBVwDdn/ZlaBcZwuBrmvDRqoOcLdbDQmo1K7/IzqQP1rElI5e3JnRncCffact0S7cWtI6oz98X7cLh1LWKmuCJRNECKD9uYYZ7XhkR6Q5EK6W+rWpFIjJFRFJEJCU7O9sDofm3Kf3akFNQzPxUPSykVnMKiuzc/cFaNqaf4q0J3RjepanRIZ3HYjbxxJD27D6azzebdJc3NcHrF7NFxAS8AfzuUmWVUlOVUolKqcSICGPPhfqCpJiGdGsZxrQVB/QvJ61GnCm2c8/MdaQePMk/xydwQ9dmRodUoRu7NKNTs1De+GG37kizBngiUWQC0eVeR7nnlQoBugDLRCQN6AV8rS9oX5qIcH+/NhzKOcP/th4xOhytljtb7GDyhymsS8vhH+MSGBHX3OiQKmUyuYZMPZRzhrkpusbtbZ5IFOuAdiISKyI2YDzwdembSqlcpVS4UipGKRUDrAZGKqVSPLDtWm9Ip0hiw+vzn5/26bs8NK8pLHEw5eMUVu0/wetj4v2i65wBHSJIbNWQfy3ZQ2GJvo7nTdVOFEopO/AwsBDYAcxVSm0TkRdEZGR111/XmU3CfX1bsyUzl1X79b3jmucV2R3c/3EqK/ce59Vb47ile5TRIV0WEVet4ujpIj5eddDocGo1j1yjUEp9p5Rqr5Rqo5R6yT3vWaXU1xWUHaBrE1fmlu4tCA+2MXW5boCneVaR3cFvPlnPT7uz+dvoroxJjL70Qj6kZ+vG9GsfwbvL9pJXWGJ0OLWWbpntBwKtZiZdG8OyXdnsOKx7z9Q8o9ju5KFPN/DjzmO8NLoL45NbGh3SVXlyaAdOnilh+soDRodSa+lE4Sfu6NWKIJuZabpWoXlAicPJI7PWs3jHUV4Y1Znbe/pvT81doxowvHNT3l9xgJO6JwOv0InCT4QF2Rif1JKvN2WRpccP1qrB7nDy29kbWLjtKM+O6MRdvWOMDqnafje0PQXFdv6jx3LxCp0o/Mg918WggBm6iq1dJbvDyeNzN/HdliP8+aZruOe6WKND8oh2kSGM7taCmb+kcfS07vzB03Si8CNRDYO4Oa4Zs9YeIvesvnCnXRmHU/H7eZv4ZlMWf7yhI5P7tjY6JI96fHB7nErx1o96yFRP04nCz0zp14aCYgefrtG3A2qXz+FUPDl/E19uzOLJYR24v38bo0PyuOhGQYxPasnstel6LBcP04nCz3RqHkrfduF88HMaRXbdyEi7NKdT8cfPN/P5+kyeGNKehwa2NTokr3lkUFssZuHNxbuNDqVW0YnCDz3Qvw3ZeUV8uSHz0oW1Os3pVDzz5RbmpmTw6PXtePT6dkaH5FVNQgOZeG0MX2zMZPfRPKPDqTV0ovBD17ZpTOfmoby3fD9O3VmgVgmlFM9+vZVZa9N5aGAbHh9cu5NEqQf6tSHYZuHvi3YZHUqtoROFHxIR7u/fhv3ZBSzeoQdv0Sr2woLtfLL6EPf3b83vh3ZARIwOqUY0rG9jct/WLNx2lE3pp4wOp1bQicJP3dilKVEN6+luPbQKbc86zQc/p3FX71Y8PbxjnUkSpe7tG0uj+jZe17UKj9CJwk9ZzCYmXxdLysGTpKTlGB2O5mPmpaZjcw/wU9eSBEBwgIUHB7RhxZ7jrNqnO9OsLp0o/NjYpGjCgqy8p2sVWjnFdidfbshkSKdIwoJsRodjmDt6taJpaCCvL9qlu+ivJp0o/FiQzcJdvWNYvOMoe4/lGx2O5iOW7DjKyTMljEn0j+7CvSXQaubR69uRevAkS3cdMzocv6YThZ+b2LsVNrOJ91foWoXmMi81g6ahgfRtp4cTHpMYRavGQby2cLe+Q7AadKLwc42DAxiTGMXn6zM5pvu4qfOOni5k2a5j3NqjBWZT3bs2cSGr+zrNjsOn+XbLYaPD8VseSRQiMlxEdonIXhF5uoL3nxCR7SKyWUSWiIj/9mnsgyZf1xq708nMX9KMDkUz2OfrM3EquK2Hfw1A5E03xzWnQ2QIb/ywG7vDaXQ4fqnaiUJEzMA7wA1AJ2CCiHS6oNgGIFEpFQfMB16t7na1c2LC6zO8S1M+Xn2Q/CK70eFoBlFKMS8lnaSYhsSG1zc6HJ9hMgm/G9qeA8cL+Gx9htHh+CVP1CiSgb1Kqf1KqWJgNjCqfAGl1FKlVGkvXauBun2VzQvu79eGvEI7s9ceMjoUzSDrD51k//ECvxvOtCYM6RRJfHQY/1y8R/eRdhU8kShaAOnlXme451XmXuB7D2xXKyc+OoxerRsxfeUBSnT1uk6al5JBkM3MTV2bGR2KzxER/jCsA1m5hfx3jf4xdaVq9GK2iNwBJAKvVfL+FBFJEZGU7OzsmgytVri/XxsO5xbyzaYso0PRatiZYjvfbMrixq7NqB9gMTocn9SnbTjXtmnMO0v3UqBP0V4RTySKTKB8XTfKPe88IjIYeAYYqZQqqmhFSqmpSqlEpVRiRIS+te9KDegQQYfIEKYu368bGNUx3285QkGxg7H6tFOVfj+sA8fzi/WNH1fIE4liHdBORGJFxAaMB74uX0BEugHv4UoSuuWLl4gIU/q1ZueRPH7arWtkdcm81HRiGgeRFNPQ6FB8WveWDRl8TST/+WkfuWf0KJGXq9qJQillBx4GFgI7gLlKqW0i8oKIjHQXew0IBuaJyEYR+bqS1WnVdHN8c5qGBvLeT7oBXl1x6MQZVu/PYUxidJ3s1+lK/W5oe/KL7Ly3fJ/RofgNj5zMVEp9B3x3wbxny00P9sR2tEuzWUzce10sL323g80Zp4iLCjM6JM3L5qemYxK4pXtV95Bopa5pFsrNcc354Oc0JvWJoUlIoNEh+TzdMrsWGp8cTUigRXcWWAc4nIr5qRlc1y6CZg3qGR2O33h8SHuKHU7eXaprFZdDJ4paKCTQyu09W/H9lsMcPFFgdDiaF/2y7zhZuYWMreMdAF6p2PD6jE2M4r9rDpFx8sylF6jjdKKope7uE4PFZOL9FQeMDkXzonkpGTSoZ2XwNZFGh+J3HhnUDgT+tWSP0aH4PJ0oaqnI0EB+1a0581LTOZFf4d3Imp/LPVPC/7Yd4VcJzQm0mo0Ox+80D6vHnb1aMT81g33Zupv+quhEUYtN6deawhInH606aHQomhd8vTmLYrtTd9lRDQ8OaEM9q5k3fthtdCg+TSeKWqxtkxAGXxPJR6vSOFus+7epbealpNOxaQidm4caHYrfahwcwL3XxfLt5sNszcw1OhyfpRNFLXd//9acPFPCvNT0SxfW/MauI3lszshlrG47UW2T+7WmQT0rf1+0y+hQfJZOFLVcYquGdG8ZxvsrDugRvmqReSnpWM3Cr7rpthPVFRpo5YH+bVi6K5uUtByjw/FJOlHUciLCxGtjOJRzhl/2nTA6HM0DShxOvtiQyfUdI2lU32Z0OLXCxGtbERESwKsLd+l+0iqgE0UdMKxzU8KCrMxep7tXrg1+3HmMEwXFjE3SbSc8Jchm4ZFBbVl7IIcVe44bHY7P0YmiDgi0mhndrQWLth0lp6DY6HC0apqXkkGTkAD6tdM9LHvS+KSWRDWsx2u6VnERnSjqiPFJLSl2OPlcDwXp147lFbJ01zFu6R6FxawPX0+yWUw8Nrg9WzJzWbjtiNHh+BT9SasjOjQNoVvLMGavS9e/lvzYlxsycTgVY3SXHV4xulsL2jYJ5vVFu3Homz/K6ERRh0xIasneY/msP3TS6FC0q6CUYm5KBt1bhtEmItjocGols0n43ZD27D2Wz5cbLhp/rc7SiaIOuSmuGfVtZmat1W0q/NHG9FPsPZavR7HzsuFdmtK1RQP+sXg3xXY9/jzoRFGn1A+wMDKhBd9uPszpQj26l7+Zl5pBoNXETXHNjA6lVhMRfje0PRknzzJH3ykI6ERR54xPiuZsiYOvN2YZHYp2Bc4WO/hmYxY3dmlGSKDV6HBqvf7tI0iOacRbP+7V3d/goRHuRGQ48E/ADLyvlHr5gvcDgI+AHsAJYJxSKs0T267tlFJQUoJyOFB2u+tRUgKl03YHyl7+tetx3usSOygniIlYEW4r3M/meWnkF3UBkwnEhJgETGYwCWIyueeXmzaZKp4vgpjNrvdFyspWNr9s2my+eL5WqYXbjpBXZNcdAF4l5XCgHA7XsWS3u16X2MF+/mtlLwH3sfZUswJeXLuFb6ad4IZrIlAOx7nP7HnHgPv4EQExVXwMVfZe2fpKj8Fyx8PVvucF1U4UImIG3gGGABnAOhH5Wim1vVyxe4GTSqm2IjIeeAUYV9V6Hbm55H71leuLUgFKgdMJyoFy2MHpcL12Olxfgk73B0FdON/pKu+eLitbbrp0viq3zLl55ZYr954q91o5XA8c56aV3T2v3PvnlXEqlN2B03H++8qp3GUUOBXihTsv7nU/p3/r8VVXixJQJgGLCWU2gdmEWNxJxWLGZCl9tmCyWDBbXc9S+rBaELMZsVpd0+75WC3ueVZ3OWu5h81VxmZDLFbEZjv3ni3APR3gmm8LAKv7fbOrW++yG8hKJ87NKPdaLpilKl6m/DF+wXqXLd5IT2sh3cmmeP/xc++ft+KLY1HKCQ7XFyHuh3K6p53uY6n0S7T0tdPpLude1v3addw4XJ9Ppx3KPt+Oc8dU2bKl5c4dY6Xvu8qrc8u4t+X60eOO0e7A6bDjdDhw2h047a5p1w8mZ9kXf/njirJphThL57mOIbmKwygI+BvAL64vtrrMEzWKZGCvUmo/gIjMBkYB5RPFKOB59/R84G0REVXFfZol6RlkPfW0B8KrHoU6d5yXHsjiPibdz04zOE3gMJ17dpjAbnY/u6ftJsFuhmLLudflyznMrrJO9zIOMzhM4ipnPrfe85Yxlf5QUeD6wYLJrFzPAmaz+1lcxybK9T3icCoEkNKdcM8v/wyuaSmdD2VlSx9SVk7KylFuedc813sCmBSYnK75olyvXdMKkwKzE8xOJxaHE4sTLA7Xw+yeNjvAYgdzMVgcCqsDLE5czw7XvNLlzE5XebPDvZ9+7H738/55hobhEU5ROMV1PLmez/982yv4jDvKjofSYwKclnPHzMVlzq3DKYLdLBetS7mPGUy4P4yuY6n0WUo/sKIQk+BqtiI4yx8rcO6YAZwXvFaq9BgSnO7NnLdMBb8rXGWk7Bg5/zgpN03F7+OFvg09kShaAOVvo8kAelZWRillF5FcoDFQaVt51bQ+hU8lUwSUmKAIRbFAsbiei3BSjKIYJ0UoisTpmqcc7vecFCo7RcpJIXaKnXaKcFCk7BQ6SyhWduzixP095tpm6Xem4P6UVC3AHIDNbCPAZCPAHIDVbMNmspZNB5ht2Ew2bBZ3OXMAVvf7NrPN9TDZCL7gtWv5i8tV9NpislxVdfPxORtZvOMoa/80mHo2Lwx6o84/EpTTiVM5cOLE6XTgdNfgHMqJE9cvSoeyY3eUUGQvpMhxliJ7IYUlZyl2FFFodz87CimyF3LWUUyRo5AiRzFF9iIKHUUUO4sodBRT7Cim0Fnses9ZTJGjhGJ7MfaSYhz2EuwlJTjtJWB3uhLPRQlJnZecKkpYpQlSwHUKAXGdXSh77XpXENeXjvuV6/9y78sF88u9X/rvWlTipNDupEE9K2aTUPaf+5/dhJRt89z6XK8dQLE4KREnJeI6XkrESYn7mCnBdbyUiAOHCM6yL2/XD5bSL3PnBa8rmy8iWMxWLBYrVrPFNW22YDNbsZitWC02AkxWbCar+3NsxWa2YjVZsYoFi8n1sJqsF01bTTYsJguB7vIWkxWryeJ+tmI127CYLVikdNr9LFYsZlvZPIt7WTGZL+s4f+jT9SzbdYzlfxhI4+AAjxwel0MphUM5yp6dyolTOc+bp1A4nO5n5cDpdPLitJYej8Uj1yg8RUSmAFMAAmMCucu0/vwC5TOwW+kXb+nDZrYRaAks+2KuZw4grHS++dz88x6WcsteZhmryerX59XHJUXzxYZMvttymFt7eKHxVul51NKXJjNmrPjSOGx2p92VVByFrmd7IUWOovMf9qJz75crV/7AVSjXs3I9Oyk3Xe790mlH6bRSrsRZbtnzyipX2VX7sqlnM9GlRWil665w2zixiKXsc2s1WwkwuaaDzDbC3J9zm8lW9gOk9MdM2XS5ZcqXKb9M+WmLyae+Ujzi8SHt+X7rYf69bB9/HtGpxrYrIljE9fe0YuwNDJ74V80Eyl9hi3LPq6hMhohYgAa4LmqfRyk1FZgK0CGug5o6ZOq5L2nTxV/YAeYATKJv3LoaPWMbERten9nrDnknUfiB0l+sQdYgo0Op1C97j/Pd4jU8Pz6BUQm6S3EjtG0SzC3do/ho9UHu7RtLswb1jA6pxnniW3Yd0E5EYkXEBowHvr6gzNfARPf0bcCPVV2fAAixhdC7eW+6R3anc+POtG3YlujQaJoENaFBQAPqWerpJFENIsK4pGjWpZ1k7zE9XrCvmpeaQUighWGdmxodSp322OB2KKX415I9RodiiGp/0yql7MDDwEJgBzBXKbVNRF4QkZHuYtOBxiKyF3gCMP4qtcat3aOwmEQ3KvJRpwtL+G7LYUYlNCfQ6ksn7eqeqIZB3N6zFXNTMtifXfd+WHnkJ7lS6julVHulVBul1Evuec8qpb52TxcqpcYopdoqpZJL75DSjBUREsDgayL5bH2m7qrABy3YdJgiu5MxPXTbCV/w8KC2BFhM/P2H3UaHUuP0uZs6bnxyNDkFxfyw/ajRoWgXmJuSTvvIYOKiGhgdigaEBwdw73WxfLv5MFszc40Op0bpRFHH9W0XQYuwenr0Ox+z91geG9NPMTYx2q/vrqtt7uvXmrAgK68t9EJjBR+mE0UdZzYJYxKjWLn3OOk5Z4wOR3Obl5KBxST8qpu+08mXhAZa+U3/Nvy0O5s1++vOGPQ6UWhl/QfNS9Hdj/uCEoeTz9ZnMrBjE8JrsIGXdnkmXhtDZGgAr9ahIVN1otBoEVaP/u0jmJuSgd2hL2ob7add2RzPL9LjTvioQKuZR69vR+rBk/y485jR4dQInSg0wNX9+JHThSzfk210KHXe3JR0woNtDOgQYXQoWiXGJkYT0ziI1xbuwlkHhkzViUID4PprIgkPtunR7wx2PL+IH3ce45buUVjN+vD0VVaziSeGdmDnkTy+3lT7x3bRn0QNcH3wb+0RxY87j3HsdKHR4dRZX27IxO5UjKmj3ar4kxFdm9GpWShv/FD7h0zViUIrMz6pJQ6nYl5qXe993xhKKealZBAfHUa7yBCjw9EuwWQSnhzWgUM5Z5hTy28E0YlCKxMbXp+esY2Ym5JeJ867+potmbnsOprH2ERdm/AXAzpEkBTTkLeW7KnVQ6bqRKGdZ0JySw6eOMPqOnSPuK+Yl5JBgMXEzfHNjQ5Fu0wiwh+Gd+RYXhEzf0kzOhyv0YlCO8/wLk0JDbQwe13trkr7msISB19tzHT//Y0de0C7MkkxjRjYIYL//LSP3LMlRofjFTpRaOcJtJq5pXsU/9t6hJMFxUaHU2cs2n6U04V23XbCTz05rCO5Z0uYunyf0aF4hU4U2kXGJUVT7HDyxYYLx5/SvGVeSjotwurRu3Vjo0PRrkKn5qGMjG/OjJVpHMurfXcN6kShXeSaZqHER4cxe92hOtNFgZEyT51l5d7j3NYjCpNJdwDor54Y0p4Sh5O3f9xrdCgepxOFVqHxSdHsPprPhvRTRodS632emoFScJtuO+HXYsLrMzYpmllrD9W6DjZ1otAqdHN8c4JsZmav1d2Pe5PT3W6ld+vGRDfy3bG7tcvz6KB2mET4Ry0b3KhaiUJEGonIDyKyx/3csIIyCSKySkS2ichmERlXnW1qNSM4wMLNcc35ZtNh8gpr550cvmBtWg6Hcs4wNknXJmqDpg0CmXRtDF9szGTXkTyjw/GY6tYongaWKKXaAUuoeCzsM8BdSqnOwHDgTREJq+Z2tRowPjmasyUOvtl02OhQaq15KRmEBFgY3rmZ0aFoHvJA/zYE2yy8vqj2DG5U3UQxCvjQPf0h8KsLCyildiul9rins4BjgO4W0w8kRIfRITKEOXr0O6/IL7Lz3ZbDjIhvTj2b2ehwNA9pWN/G/f1b88P2o6w/dNLocDyiuokiUilV+nPzCBBZVWERSQZsQIU3G4vIFBFJEZGU7Gzd3bXRRITxydFsyshle9Zpo8Opdb7dnMXZEgdjdJcdtc7dfWIJD7bx6v921oo7By+ZKERksYhsreAxqnw55fprVPoXEZFmwMfA3UqpCrtaVEpNVUolKqUSIyJ0pcMXjO7WApvFpGsVXjA3JYM2EfXpFh1mdCiah9UPsPDwwLas3p/Dij3HjQ6n2i6ZKJRSg5VSXSp4fAUcdSeA0kRQ4XBPIhIKfAs8o5Ra7ckd0LwrLMjGDV2a8sWGTApLam+nZzVtX3Y+qQdPMjYxGhHddqI2mtCzJS3C6vFaLRgytbqnnr4GJrqnJwJfXVhARGzAF8BHSqn51dyeZoBxSdGcLrTz/VZ9UdtT5qdmYDYJo7u3MDoUzUsCLGYeH9KeLZm5fL/1iNHhVEt1E8XLwBAR2QMMdr9GRBJF5H13mbFAP2CSiGx0PxKquV2tBvVu3ZiYxkF69DsPsTucfJaawYD2ETQJCTQ6HM2LRndrQbsmwby+aJdfj0dfrUShlDqhlLpeKdXOfYoqxz0/RSk12T39iVLKqpRKKPfY6IHYtRoiIoxNimbtgRz2Z+cbHY7fW7HnOMfyihijOwCs9cwm4XdDO7A/u4DP1/tv32m6ZbZ2WW7rEYXZJMzR3Y9X29yUdBrVtzGoYxOjQ9FqwLDOkcRHh/Hm4t1+e51PJwrtsjQJCeT6jk34bH1GrR8f2JtyCopZvONo2d1kWu0nIjw1rANZuYV8svqg0eFcFf1J1S7bhOSWHM8vZsmOo0aH4re+2phJiUPpthN1zLVtw7mubTjvLttHfpHd6HCumE4U2mXr1z6CZg0C9eh31TA3JYOuLRrQsWmo0aFoNezJYR3IKSjm/RX7jQ7lilmMDkDzH2aTMCYxmrd+3EPGyTNENdS9nV6JrZm57Dh8mhdHdfbqdkpKSsjIyKCwsPYNoOPPbMAnt0ZRZM9n27bt1R57JDAwkKioKKxW7w+dqxOFdkXGJkbx1o97mJeSweND2hsdjl+Zn5qBzWJiZLx3205kZGQQEhJCTEyMbsznYwpLHOw5mkfD4ACah9W76vUopThx4gQZGRnExsZ6MMKK6VNP2hWJahhE33YRzEtJx+H079amNanI7uDLjZkM7RRJgyDv/gIsLCykcePGOkn4oECrmbAgGycKiqt1U4iI0Lhx4xqrNepEoV2x8UnRZOUWsnyP7rjxci3efoxTZ0oYW0NtJ3SS8F2RoQEA1R5buyb/jXWi0K7Y4GsiaVzfpke/uwJzU9Jp1iCQPm3DjQ5FM5jNYqZxfRsnC4r9pl2FThTaFbNZTNzaI4olO45V+1dRXXAkt5AVe7LLGi3WBS+99BKdO3cmLi6OhIQE1qxZc8XrmDBhAnFxcfzjH//wQoSVmzRpEvPne7ZbupkzZ5KVlVX2+s+/e5j9e3Zx9LR/HD/6YrZ2VcYlRTN1+X4+S83kNwPaGB2OT/tsfQZO5WrdXhesWrWKBQsWsH79egICAjh+/DjFxcWXvbzdbuf48eOsW7eOvXv3ejHSmjNz5ky6dOlC8+bNAZgxfTpHcgs5llfI2WI79Wy+/VWsaxTaVWkTEUxyTCPmrDvk910oe5NSinkp6STHNqJV4/pGh1MjDh8+THh4OAEBrnPx4eHhZV+QMTExHD/uGp8hJSWFAQMGAPD8889z55130qdPH+68806GDh1KZmYmCQkJrFixgmnTppGUlER8fDy33norZ86cAeDo0aOMHj2a+Ph44uPj+eWXXwD45JNPSE5OJiEhgfvvvx+H4+JTPKmpqfTv358ePXowbNgwDh++uHfkisrs3LmT5OTksjJpaWl07doVgBdeeIGkpCS6dOnClClTUEoxf/58UlJSuP3220lISODs2bMMGDCAg7s2YzYJU2d+TNeuXenSpQtPPfVU2XqDg4N55plniI+Pp1evXhw9alxDV99OY5pPG58czRNzN7F6fw692zQ2OhyflHLwJGknzvDwoHaGbP//vtnm8dEJOzUP5bmbK28LMnToUF544QXat2/P4MGDGTduHP3797/kerdv387KlSupV68eaWlpjBgxgo0bN7q22akT9913HwB//vOfmT59Oo888giPPvoo/fv354svvsDhcJCfn8+OHTuYM2cOP//8M1arlQcffJBPP/2Uu+66q2xbJSUlPPLII3z11VdEREQwZ84cnnnmGWbMmHFZZYqLizlw4ACxsbHMmTOHcePGAfDwww/z7LPPAnDnnXeyYMECbrvtNt5++21ef/11EhMTy9ZvNplw5p/g5f/7f/yyZh3RTSMYOnQoX375Jb/61a8oKCigV69evPTSS/zhD39g2rRp/PnPf778fygP0jUK7ard0KUZIYEWPfpdFeauS6e+zcyNXZsaHUqNCQ4OJjU1lalTpxIREcG4ceOYOXPmJZcbOXIk9epV3LZg69at9O3bl65du/Lpp5+ybds2AH788Ud+85vfAGA2m2nQoAFLliwhNTWVpKQkEhISWLJkCfv3n98aeteuXWzdupUhQ4aQkJDAX/7yFzIyMi67zNixY5kzZw7AeYli6dKl9OzZk65du/Ljjz+WxVmZvds2k3xtXxy2EMxmM7fffjvLly8HwGazMWLECAB69OhBWlraJf+G3qJrFNpVq2czM7pbC2avS+f5M8WEBdmMDsmnFBTZ+XbLYW6Oa06QQeegq/rl701ms5kBAwYwYMAAunbtyocffsikSZOwWCw4na72Axe2Aahfv/JTc5MmTeLLL78kPj6emTNnsmzZskrLKqWYOHEif/vb36os07lzZ1atWnVVZcaNG8eYMWO45ZZbEBHatWtHYWEhDz74ICkpKURHR/P8889fsp2DySQEWk2cKbaTV3h+H1BWq7XsFliz2YzdblwfUbpGoVXLuKRoiu1Ovtzgv33te8t3Ww5zpthR5zoA3LVrF3v27Cl7vXHjRlq1agW4rlGkpqYC8Nlnn132OvPy8mjWrBklJSV8+umnZfOvv/56/v3vfwPgcDjIzc3l+uuvZ/78+Rw75hqZOScnh4MHz++1tUOHDmRnZ5clgZKSkot+/VdVpk2bNpjNZl588cWy2kRpUggPDyc/P/+8O6dCQkLIy8u7aL+Sk5NZ/fNKzuSeIvNkAbNmzbqs03Q1TScKrVo6N29AXFQDZq9L1xe1LzAvJYPW4fXp0aqh0aHUqPz8fCZOnEinTp2Ii4tj+/btPP/88wA899xz/Pa3vyUxMRGz2XzZ63zxxRfp2bMnffr0oWPHjmXz//nPf7J06VK6du1Kjx492L59O506deIvf/kLQ4cOJS4ujiFDhlx0odpmszF//nyeeuop4uPjSUhIKLsQfrllxo0bxyeffMLYsWMBCAsL47777qNLly4MGzaMpKSksrKTJk3igQceKLuYXapZs2a8/PLL3DN2BDcP7E3nuARGjRp12X+XmiLVObhFpBEwB4gB0oCxSqmTlZQNBbYDXyqlHr7UuhMTE1VKSspVx6bVnE/XHOSZL7by5UN9SIgOMzocn5B2vIABry/jD8M78OCAtjW67R07dnDNNdfU6Da16lFKsfdYPg6laB8ZgukyW11X9G8tIqlKqcRKFrkq1a1RPA0sUUq1A5a4X1fmRWB5Nben+aCR8c2pZzXrltrlzE/NwCRwa/e6ddpJuzoiQmSDQIrtTk4WXH6bk5pS3UQxCvjQPf0h8KuKColIDyASWFTN7Wk+KCTQyoi4Zny9KcsvB2XxNIdTMT81g37tI4gMDTQ6HM1PhARYqG+zcDSvyOc63KxuoohUSpWe/DuCKxmcR0RMwN+B319qZSIyRURSRCQlO1t3OOdPxie35EyxgwWbsi5duJZbufc4R04X1lgHgFrtICI0bRCI3eHkREGR0eGc55KJQkQWi8jWCh7nXXFRrosdFaXBB4HvlFIZFbx3HqXUVKVUolIqMSIi4rJ3QjNe95ZhtGsSrEe/w9UBYFiQleuvaWJ0KJqfqR9gITTQSnZeEXaH74xNf8mbu5VSgyt7T0SOikgzpdRhEWkGHKugWG+gr4g8CAQDNhHJV0pVdT1D8zMiwvjklry4YDs7j5yus0N9njpTzA/bjvLrni0JsFz+XT2aVioyNJA9x/LIzi+iWYOrH9zIk6p76ulrYKJ7eiLw1YUFlFK3K6VaKqVicJ1++kgnidppdLcW2MwmZq+tu7WKrzdlUexw1rm2E5rn1LO5BzfKL6bER2oV1U0ULwNDRGQPMNj9GhFJFJH3qxuc5l8a1bcxrEtTvtiQ6Tf97Hva3JR0OjULpXPzBkaHYpi0tDS6dOly3rznn3+e119/vdrrvtr1ZGVlcdtttwGuBoDfffddtWPxpsjQAJSCY6d941pFtRKFUuqEUup6pVQ7pdRgpVSOe36KUmpyBeVnXk4bCs1/jU+KJvdsCQu3HTE6lBq3Pes0WzNPM1bXJnxO8+bNy1pK+0OiCLCYaVTfSk5BMUV243906ZbZmkf1bt2Ylo2CmFUH21TMS03HZjYxKqGF0aH4hdzcXFq1alXW91NBQQHR0dGUlJSwb98+hg8fTo8ePejbty87d+68aPmNGzfSq1cv4uLiGD16NCdPutr67t27l8GDBxMfH0/37t3Zt29fWS2nuLiYZ599ljlz5pCQkMCcOXNo164dpXdZOp1O2rZtiy/cddkkNBAR36hV6E4BNY8ymYRxSdG8tnAXB44XEBteN8ZgKLY7+WpjFoM7NaFhfR/qHPH7p+HIFs+us2lXuOHlaq+mQYMGJCQk8NNPPzFw4EAWLFjAsGHDsFqtTJkyhf/85z+0a9eONWvW8OCDD/Ljjz+et/xdd93FW2+9Rf/+/Xn22Wf5v//7P958801uv/12nn76aUaPHk1hYSFOp7Os3yebzcYLL7xASkoKb7/9NgA7d+7k008/5bHHHmPx4sXEx8fjC3ddWs0mGgfbyM4rIjwkgHpW426O0DUKzeNKh/ycU4dulf1x51FyCooZo9tOlPV4ejnzx40bV9Zd9+zZsxk3bhz5+fn88ssvjBkzpmzgoQv7asrNzeXUqVNlHehNnDiR5cuXk5eXR2ZmJqNHjwYgMDCQoKCgKuO95557+OijjwCYMWMGd99995XtsBdFBAdgNglHc40dMlXXKDSPiwwNZGCHJsxPzeB3Q9tjNdf+3yNzUzKIDA2gXzvjf4mexwO//K9U48aNy04DlcrJySE2NvaisiNHjuRPf/oTOTk5pKamMmjQIAoKCggLCysbtMjboqOjiYyM5Mcff2Tt2rXn9U5rNIvZRERwAEdOF1JQZKd+gDFf2bX/CNYMMT4pmuP5RSzZUVHTmtrl2OlClu06xq3dXTWpui44OJhmzZqVnSrKycnhf//7H9ddd12FZZOSkvjtb3/LiBEjMJvNhIaGEhsby7x58wBXh3mbNm06b7kGDRrQsGFDVqxYAcDHH39M//79CQkJISoqii+//BKAoqKismFTS1XU5ffkyZO54447GDNmzBX1alsTGgcHYDGZOHK60LAemnWi0LxiQIcIIkMD6sTod59vyMSpXKfcNJePPvqIF198kYSEBAYNGsRzzz1HmzZtKixb2l136bgOAJ9++inTp08nPj6ezp0789VXFzXR4sMPP+TJJ58kLi6OjRs3lg1B+vHHH/Ovf/2LuLg4rr32Wo4cOf8OvIEDB7J9+/ayi9ngqtnk5+f71GmnUmaT0CQ0gIIiu2F9qVWrm3Fv0t2M+7/XF+7i3WV7WfnUIJqH+UYLU09TSnH9Gz/RKMjG/N9ca3Q4gO5m/GqkpKTw+OOPl9VQfI1TKXYfycNsEto2CS673uMv3YxrWqXGJUXjVK4BfGqr9YdOsT+7QHcA6Mdefvllbr311iqHTjWaSYTI0EDOljjIPVtS89uv8S1qdUZ0oyCuaxvO3JR0n+s22VPmpaRTz2rmxrhmRoeiXaWnn36agwcPVngNxZeEBVkJtJo5erqoxq9V6EShedX45GgyT51l5d7jRoficWeK7SzYfJib4poRbNDdKFrdIe5aRZHdwckzNTu4kU4UmlcN6RRJwyBrrRz97n9bj5BfZGeMvoit1ZDQQAtBNgtHTxfhrMFauk4UmlcFWMzc2j2KH7Yf5Xi+8V0ReNLclHRaNQ4iObaR0aFodYSI0DQ0gBKHkxM1OGSqThSa141PjsbuVHyWWnsuah86cYbV+3MY0yOq0pbImuYNwYFWggMsZOcV4qyhaxU6UWhe17ZJCImtGjJnXbphDYY8bX5qOiJwS3d92ulCAwcOZOHChefNe/PNN/nNb35T7XUHBwdf1XLPPvssixcvLovlwkZ4/qZpg0DsTkV+Yc20q9CJQqsR45Ki2X+8gLUHcowOpdqcTsVn6zO5rm14rW0fUh0TJkxg9uzZ582bPXs2EyZMMCgieOGFFxg82DVYZ21IFEE2Cw3qWckvsnOiBk7p6kSh1Yib4poREmCpFR0F/rLvBJmnzuq2E5W47bbb+Pbbbykudp1DT0tLIysri759+55X7umnn+add94pe11+UKLXXnuNpKQk4uLieO655y7ahlKKJ598ki5dutC1a9eyFtYAr7zyCl27diU+Pp6nn3YNpjlp0iTmz5/Pv/71L7Kyshg4cCADBw5kxowZPPbYY2XLTps2jccff9xjfwtvigwNRCl4d9k+r2+rWvf0iUgjYA4QA6QBY5VSJyso1xJ4H4gGFHCjUiqtOtvW/EuQzcLIhObMT83guZs70yDIanRIV6XE4eSjVWmEBloY0inS6HAu6ZW1r7Az5+KxHKqjY6OOPJX8VKXvN2rUiOTkZL7//ntGjRrF7NmzGTt27EXXcsaNG8djjz3GQw89BMDcuXNZuHAhixYtYs+ePaxduxalFCNHjmT58uX069evbNnPP/+cjRs3smnTJo4fP05SUhL9+vVj48aNfPXVV6xZs4agoCBycs6vwT766KO88cYbLF26lPDwcPLz83nppZd47bXXsFqtfPDBB7z33nse/Gt5T6DVTFCAmY9XH+Se62Jp4cXabXVrFE8DS5RS7YAl7tcV+Qh4TSl1DZAM1P6e4rSLTEhuSZHdyVebMo0O5YrlFZYwbfl++r26lEXbj/Lrnq0INHB8AF9X/vRTZaedunXrxrFjx8jKymLTpk00bNiQ6OhoFi1axKJFi+jWrRvdu3dn586d7Nmz57xlV65cyYQJEzCbzURGRtK/f3/WrVvH4sWLufvuu8u6Fm/UqOo70oKDgxk0aBALFixg586dlJSU0LVrVw/9FbwvJNAKCv65eLdXt1PdVkKjgAHu6Q+BZcB5PzVEpBNgUUr9AKCUyq/mNjU/1aVFAzo3D2XW2nTu7NXKL+4WOpx7lpk/p/HfNYfIK7LTM7YRL43uwoD2TYwO7bJU9cvfm0aNGsXjjz/O+vXrOXPmDD169Kiw3JgxY5g/fz5Hjhwp6xRQKcUf//hH7r///hqJdfLkyfz1r3+lY8eOPtkpYFUsJuGOXq2Y+csBpvRrQ9smV3ex/1KqW6OIVEqVjihyBKioLt4eOCUin4vIBhF5TUT0T7E6anxyS3YcPs2WzFyjQ6nSjsOneWLuRvq+spRpK/bTv0MEXz3Uhzn392ZQx0hMujvxKgUHBzNw4EDuueeeKi9ijxs3jtmzZzN//nzGjBkDwLBhw5gxYwb5+a7flJmZmWUj1JXq27cvc+bMweFwkJ2dzfLly0lOTmbIkCF88MEHZRerLzz1BBd3M96zZ0/S09P573//a+gF96v10MA21LOaeeOHXV7bxiVrFCKyGGhawVvPlH+hlFIiUtG9jxagL9ANOITrmsYkYHoF25oCTAFo2bLlpULT/NCohOa89O12Zq1NJy4qzOhwzqOUYuXe40xdvp8Ve44TZDNzZ+9W3NMnluhGVY+Spl1swoQJjB49+qI7oMrr3LkzeXl5tGjRgmbNXP1lDR06lB07dtC7d2/AlXQ++eQTmjQ5V4sbPXo0q1atIj4+HhHh1VdfpWnTpgwfPpyNGzeSmJiIzWbjxhtv5K9//et525wyZQrDhw+nefPmLF26FICxY8eyceNGGjZs6Ok/g9c1Dg7g3r6t+deSPWzJ8M4PsGp1My4iu4ABSqnDItIMWKaU6nBBmV7AK0qp/u7XdwK9lFIPVbVu3c147fW7uZv439bDrH1msGEjdpVX4nCyYHMWU5cfYMfh00SEBDDp2hhu79mSsCAfGv/6Muluxq/ciBEjePzxx7n++uuNDuWKlP5b5xWW0O/VpXRp0YBPJvfyuW7GvwYmuqcnAhePLgLrgDARKR0jchCwvZrb1fzY+ORoCoodfLv58KULe1H5C9SPz9mE3eHk1dviWPnUQB4a2NYvk4R2ZU6dOkX79u2pV6+e3yWJ8kICrTw4oC0r9nin883q/px7GZgrIvcCB4GxACKSCDyglJqslHKIyO+BJeK6epkKTKvmdjU/ltiqIW0i6jN73SHGJtV8W4SsU2eZ+Usas9wXqHu3bsxfR3elf/sIfe2hjgkLC2P3bu/eMVRT7uzdihk/H+CgF9ZdrUShlDoBXJSGlVIpwORyr38A4qqzLa32EBHGJ7Xkpe92sPtoHu0jQ2pku9uzTjNtxX6+2ZTlaszTtRlT+rama1SDGtm+pnlToNXMb69vx2ovrNv4E8RanXRL9xa8unAns9em8+zNnby2HaUUK/YcZ9qKcxeo7+odw919YvQFaq3WGZMYjTfu29KJQjNE4+AAhnZqyucbMvjD8A4eb7xWbC+9QL2fnUfyaBISwFPDO/Lr5JZ+2ypc0y7F7KVTpzpRaIYZnxzNt1sOs2j7UUbGN/fIOk8XljBrzSE++DmNI6cLaR8ZzGu3xTEyoTkBFt18R9Ouhu4UUDNMnzbhRDWs55HR77JOneWlb7dz7d9+5G/f76R1RH0+uDuJhY/1Y0xitE4SNezIkSOMHz+eNm3a0KNHD2688cYav2i8bNkyfvnllxrdZkXS0tLo0qWL0WFUi65RaIYxmYRxidH8/YfdHDxRQKvG9a94Hduycpm2fD8LNh9GASPimnFf39Z0aaEvUBtFKcXo0aOZOHFiWWO7TZs2cfToUdq3b39Z63A4HJjN5kpfX45ly5YRHBzMtddee0XLaRfTNQrNULclRmESrqj7caUUP+3O5o7313DTv1byw/ajTLw2huV/GMg/x3fTScJgS5cuxWq18sADD5TNi4+Pp2/fvixbtowRI0aUzX/44YeZOXMmADExMTz11FN0796defPmXfR60aJF9O7dm+7duzNmzJiyLj5iYmJ47rnn6N69O127dmXnzp2kpaXxn//8h3/84x8kJCSwYsWK82L86aefSEhIICEhgW7dupGXl8eyZcvo168fN910Ex06dOCBBx7A6XQCVLrt1NRU+vfvT48ePRg2bBiHDx8umx8fH098fPx5Xan7K12j0AzVrEE9BnZowrzUDJ4Y0h6LufLfLsV2J99symLaCtcF6sjQAJ6+oSMTklvSoJ6+QF2RI3/9K0U7PNvNeMA1HWn6pz9V+v7WrVsr7QTwUho3bsz69esB13gVpa+PHz/OLbfcwuLFi6lfvz6vvPIKb7zxBs8++ywA4eHhrF+/nnfffZfXX3+d999/nwceeIDg4GB+//vfX7Sd119/nXfeeYc+ffqQn59PYGAgAGvXrmX79u20atWK4cOH8/nnnzNgwAD+8pe/XLTtP/7xjzzyyCN89dVXREREMGfOHJ555hlmzJjB3Xffzdtvv02/fv148sknr+pv4Ut0otAMNy4pmiU7j7F0V3aFYzzkni1h1tpDfPDzAY6eLqJDZAivj4lnZHxzbBZdKa5NSnuQvfD16tWr2b59O3369AGguLi4rC8ogFtuuQWAHj168Pnnn19yO3369OGJJ57g9ttv55ZbbiEqyjWkbXJyMq1btwZcfVWtXLmSwMDACre9a9cutm7dypAhQwDX6bFmzZpx6tQpTp06VTZ+xp133sn3339/1X8TX6AThWa4QR2b0CQkgNlrD52XKDJPneWDlQeYvS6d/CI7fdo25tXb4unXLtwvuij3BVX98veWzp07M3/+/Arfs1gsZadzAAoLC897v379+hW+VkoxZMgQZs2aVeF6AwICADCbzdjtlx5H+umnn+amm27iu+++o0+fPmVjfF/4uRKRSre9ZcsWOnfuzKpVq86bf+rUqUtu39/on2Oa4SxmE7f1iGLprmMcyS1ka2Yuv529gX6vLuWDX9IYfE0TFjxyHZ9O7kX/9hE6Sfi4QYMGUVRUxNSpU8vmbd68mRUrVtCqVSu2b99OUVERp06dYsmSJZe1zl69evHzzz+zd+9eAAoKCi55F9WF3YmXt2/fPrp27cpTTz1FUlISO3e6Ts+tXbuWAwcO4HQ6mTNnDtddd12l2+7QoQPZ2dlliaKkpIRt27YRFhZGWFgYK1euBODTTz+9rH30ZTpRaD5hXFI0TgW3/vsXRry1kiU7jnFPH9cF6jf1BWq/IiJ88cUXLF68mDZt2tC5c2f++Mc/0rRpU6Kjoxk7dixdunRh7NixdOvW7bLWGRERwcyZM5kwYQJxcXH07t277Mu9MjfffDNffPFFhRez33zzTbp06UJcXBxWq5UbbrgBgKSkJB5++GGuueYaYmNjGT16dKXbttlszJ8/n6eeeor4+HgSEhLKbsf94IMPeOihh0hISKA6PXT7imp1M+5Nupvxuue+j1LYkpHL3X1imNCzJaGB+gL11dDdjF+dZcuW8frrr7NgwQKjQ7lsFf1bi4jHuxnX1yg0nzH1TtedMvrUkqb5Fp0oNJ+hE4RmpAEDBjBgwACjw/BJ+hqFptVCvnpKWfOcmvw31olC02qZwMBATpw4oZNFLaaU4sSJE2UNBb1Nn3rStFomKiqKjIwMsrOzjQ5F86LAwMCyhoLeVq1EISKNgDlADJAGjFVKnayg3KvATbhqMD8Av1X6546meYXVaiU2NtboMLRapLqnnp4Gliil2gFL3K/PIyLXAn1wDYXaBUgC+ldzu5qmaVoNqW6iGAV86J7+EPhVBWUUEAjYgADAChyt5nY1TdO0GlLdRBGplDrsnj4CXNSjm1JqFbAUOOx+LFRK7ajmdjVN07QacslrFCKyGGhawVvPlH+hlFIictF1BxFpC1wDlF51+UFE+iqlVlRQdgowxf2ySES2Xiq+GtAAyDU6CA/Q++F7asu+1Jb9CAeOGx2EB3Tw9AovmSiUUoMre09EjopIM6XUYRFpBhyroNhoYLVSKt+9zPdAb+CiRKGUmgpMdZdL8XQz9KshIlOVUlMuXdK36f3wPbVlX2rRfvjEd051iYjH+z6q7qmnr4GJ7umJwFcVlDkE9BcRi4hYcV3I9qdTT98YHYCH6P3wPbVlX2rLfmiVqFangCLSGJgLtAQO4ro9NkdEEoEHlFKTRcQMvAv0w3Vh+39KqScuY921IrtrmuYfast3jjf2w2d7jxWRKe5TUZqmaV5XW75zvLEfPpsoNE3TNN9Qp/p6EpHhIrJLRPaKyNPueYNEZL2IbBWRD0Wkwgv8IjJRRPa4HxPLze8hIlvc6/yX1EAXqCIyQ0SOlb8rTEReE5GdIrJZRL4QkbBKlr3ob+CeHysia9zz54iIzaD9eF5EMkVko/txo6/vRxX7kiAiq937kSIiyZUs6xOfLRGJFpGlIrJdRLaJyG/d88e4Xzvdp5UrW96n/k00D1JK1YkHYAb2Aa1xNf7bBHQC0oH27jIvAPdWsGwjYL/7uaF7uqH7vbVAL0CA74EbamBf+gHdga3l5g0FLO7pV4BXLvdv4H5vLjDePf0f4DcG7cfzwO+v5t/SqP2oYl8WlX4egBuBZb782QKaAd3d0yHAbvcxcg2uWy6XAYn+8m+iH5571KUaRTKwVym1XylVDMwGbgWKlVKlg+/+4J53oWHAD0qpHOXqy+oHYLj7luBQpdRq5ToKPqLi1ukepZRaDuRcMG+RUqp0VPnVnGu3Ul5Ff4NR7l+qg4D57nKVtbL3qIr24zL51H5ApfuigFD3dAMgq4JFfeazpZQ6rJRa757Ow3V3Ygul1A6l1K5LLO5z/yaVnEG4rNqNiPzRXWaXiAyrap0G7cd0EdnkPoMwX0SCvbkfdSlRtMBVeyiVgashoaVcdfo2IBpARBJF5P0qlm3hfmRUMN9o9+D6BYqINBeR79zzK9uPxsCpconG6P142H0AzBCRhuC3+/EY8JqIpAOvA38E//hsiUgM0A1YU0UZn/03Edfdlu8AN+CqFU0QkU64atv/UEq1BU4C91awbCdgPNAZGA68KyLmKtZpxH48rpSKV0rF4WqC8LA396MuJYqKKFx/yH+IyFogD3AAKKVSlFKTjQzuaojIM4Ad+BRAKZWllKrwPL+P+jfQBkjA1eXL38Ev9wPgN7gO6GjgcWA6+P5ny/3r9DPgMaXU6crK+fi/SYU1HC6vdjMKmK2UKlJKHQD2utdX2Tq9qcJtlv67uGts9XB9l3ltP+pSosjEXVtwiwIylVKrlFJ9lVLJwHJc52Uva1n3I6qC+YYQkUnACOB29+mKC1W2HyeAMDl3Id+w/VBKHVVKOZRSTmAarg/1hXx+P9wmAp+7p+dxZftiyGdLXI1iPwM+VUp9fqny5fjav0llNZwKazciMlJEXrjEspXN96ZKtykiH+DqY68j8JZ7nlf2oy4linVAO/c5ShuumsTXItIEQEQCgKdwXWy70EJgqIg0dJ8KGYqrc8PDwGkR6eXO7HdRcet0rxOR4cAfgJFKqTOVFKvwb+BOKktxnXqDylvZe5373Hyp0UBF/X35/H64ZXGuS/1BwJ4KyvjMZ8u9nenADqXUG1e4uL/8m1RIKfW1UupZo+O4Ekqpu4HmuK4ljXPP88p+1JlE4f4V8TCuA3MHMFcptQ14UkR2AJuBb5RSP8L555GVUjnAi7gOhnXAC+55AA8C7+Oq1u3DfW3Am0RkFrAK6CAiGSJyL/A2rjtVfhDX7Zj/cZctO49cxd8AXEnyCRHZi+u88nSD9uNVcd0SuhkYiOuUjU/vRxX7ch/wdxHZBPwVd4eXPvzZ6gPcCQyScrcni8hoEcnA1UfbtyKy0L0fvvxvUlkN53JqN1XV8iqa701VblMp5eDcjTmXu+wV74ducKdpWq3jTga7getxfQmuA36N6/brz5RSs90/pjYrpd69YNnOwH9xnSpsjmtQtna4blO+aJ3lEmJN7keRUmqvuxb4GoBS6vfe2g89ZramabWOUsouIqU1HDMwQym1TUSeAmaLyF+ADbhrNyIyElcbkWfd5eYC23HdGPKQ+5c7Fa2zpvcDV41thYiE4vrS34Trxgmv7YeuUWiapmlVqjPXKDRN07SroxOFpmmaViWdKDRNqxWk4q4uVpS7gytLRL6sYvlfiYgSkY41FrSf0IlC0zS/V1m3FO7GtAlKqQRcty9X1YhwArDS/ayVoxOFpmm1QZXdUrjvEBoEfFnRwu5uS67D1ffT+HLzB4jIgnKv33b3gIC7jclOEUkVVzfwCy5cb22hE4WmabXBpbql+BWwpIq+q0bhGqZ5N3BCRHpUtTERCQTew9X1ew8g4moD9wc6UWiaVhdMAGZd4v3Z7unZXPr0U0dgv7uzPS6xbr+nG9xpmlYbVNothYiE4zo1Nbr0TXc3JJFACq4+0gYBXUVE4WqEpkTkSVwN1cr/oA704j74LF2j0DStNqiwU0L3e7cBC5RShaWFlVLD3Be5J7vf/1gp1UopFePuFv4A0Bc4CHQSkQBxDS98vXsVu4DW4hq3A9yd8tVWOlFomub3LtEp4XgufdrpiwvmfQZMUEql4xrKdav7eYN7e2dxddr4PxFJxTWWTa5n9sb36C48NE3TroKIBCul8t0d870D7FFK/cPouLxB1yg0TdOuzn0ishHYhmtM9PeMDcd7dI1C0zRNq5KuUWiapmlV0olC0zRNq5JOFJqmaVqVdKLQNE3TqqQThaZpmlYlnSg0TdO0Kv1/fJZkxIPTowgAAAAASUVORK5CYII=",
      "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, 4493.36it/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\\jan\\code\\mikeio\\mikeio\\dfsu.py:2709: FutureWarning: setting items is deprecated, please supply data in the form of a Dataset\n",
      "  warnings.warn(\n",
      "c:\\users\\jan\\code\\mikeio\\mikeio\\dfsu.py:2715: FutureWarning: supplying data as a list of numpy arrays is deprecated, please supply data in the form of a Dataset\n",
      "  warnings.warn(\n",
      "c:\\users\\jan\\code\\mikeio\\mikeio\\dfsu.py:2752: 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, 1005.59it/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",
    "* Renamed variables\n",
    "\n",
    "First inspect the source file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 9/9 [00:00<00:00, 1827.94it/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": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ds = dfs.read()\n",
    "ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "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:  eastward_sea_water_velocity <u velocity component> (meter per sec)\n",
       "  1:  northward_sea_water_velocity <v velocity component> (meter per sec)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from mikeio.eum import ItemInfo, EUMType\n",
    "from mikeio import Dataset\n",
    "\n",
    "sourcefilename = filename\n",
    "outfilename2 = \"HD2D_selected.dfsu\"\n",
    "\n",
    "newds = ds[[\"U velocity\", \"V velocity\"]].rename({'U velocity': 'eastward_sea_water_velocity',\n",
    "                                                 'V velocity': 'northward_sea_water_velocity'})\n",
    "newds\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 9/9 [00:00<00:00, 9413.65it/s]\n"
     ]
    }
   ],
   "source": [
    "dfs.write(outfilename2, newds)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Read the newly created file to verify the contents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 9/9 [00:00<00:00, 8964.32it/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:  eastward_sea_water_velocity <u velocity component> (meter per sec)\n",
       "  1:  northward_sea_water_velocity <v velocity component> (meter per sec)"
      ]
     },
     "execution_count": 24,
     "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": 25,
   "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": 26,
   "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.9.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}