{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "28db61ce-3915-403a-9d2d-373e6a484d7e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import required libraries\n",
    "import rasterio\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "32352e5f-3cba-4366-bba2-cc05a0f0ed7a",
   "metadata": {},
   "source": [
    "<hr/>\n",
    "\n",
    "### Open raster file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "1bf07874-a7de-48f8-b92d-6a4eeca1317f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lakeRst = rasterio.open('../Rst/lakeBottomElevation.tif')\n",
    "lakeRst.count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2a030320-1eef-4833-9f4f-6fe5ccbb6318",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.5, 0.5)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#raster resolution\n",
    "lakeRst.res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "011ed07a-bb42-48bb-a055-a973a7132c9c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'PROJCS[\"NAD83(2011) / Puerto Rico and Virgin Is.\",GEOGCS[\"NAD83(2011)\",DATUM[\"NAD83_National_Spatial_Reference_System_2011\",SPHEROID[\"GRS 1980\",6378137,298.257222101,AUTHORITY[\"EPSG\",\"7019\"]],AUTHORITY[\"EPSG\",\"1116\"]],PRIMEM[\"Greenwich\",0,AUTHORITY[\"EPSG\",\"8901\"]],UNIT[\"degree\",0.0174532925199433,AUTHORITY[\"EPSG\",\"9122\"]],AUTHORITY[\"EPSG\",\"6318\"]],PROJECTION[\"Lambert_Conformal_Conic_2SP\"],PARAMETER[\"latitude_of_origin\",17.8333333333333],PARAMETER[\"central_meridian\",-66.4333333333333],PARAMETER[\"standard_parallel_1\",18.4333333333333],PARAMETER[\"standard_parallel_2\",18.0333333333333],PARAMETER[\"false_easting\",200000],PARAMETER[\"false_northing\",200000],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]],AXIS[\"Easting\",EAST],AXIS[\"Northing\",NORTH],AUTHORITY[\"EPSG\",\"6566\"]]'"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#all raster crs info\n",
    "lakeRst.crs.wkt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "913ddbd3-ca99-4d22-a768-c77e614da5f6",
   "metadata": {},
   "source": [
    "<hr/>\n",
    "\n",
    "### Get raster values as Numpy array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "de1e5743-8edf-4850-ba28-97287de132d3",
   "metadata": {},
   "outputs": [],
   "source": [
    "lakeBottom = lakeRst.read(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "72435881-649a-4794-8cbe-0ae8e0171d5c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-3.4028235e+38, -3.4028235e+38, -3.4028235e+38, -3.4028235e+38,\n",
       "        -3.4028235e+38],\n",
       "       [-3.4028235e+38, -3.4028235e+38, -3.4028235e+38, -3.4028235e+38,\n",
       "        -3.4028235e+38],\n",
       "       [-3.4028235e+38, -3.4028235e+38, -3.4028235e+38, -3.4028235e+38,\n",
       "        -3.4028235e+38],\n",
       "       [-3.4028235e+38, -3.4028235e+38, -3.4028235e+38, -3.4028235e+38,\n",
       "        -3.4028235e+38],\n",
       "       [-3.4028235e+38, -3.4028235e+38, -3.4028235e+38, -3.4028235e+38,\n",
       "        -3.4028235e+38]], dtype=float32)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "#raster sample\n",
    "lakeBottom[:5,:5]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c0ff0e20-3e8a-4cf3-8b80-fc3df768d46e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 122688 characters]\n",
      "text/plain": [
       "<Figure size 864x864 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#replace value for np.nan\n",
    "noDataValue = np.copy(lakeBottom[0,0])\n",
    "\n",
    "lakeBottom[lakeBottom==noDataValue]= np.nan\n",
    "plt.figure(figsize=(12,12))\n",
    "plt.imshow(lakeBottom)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e6196b0d-ce1f-4ee2-866d-875d81cdf4b8",
   "metadata": {},
   "source": [
    "<hr/>\n",
    "\n",
    "### Lake volume calculation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7e377e92-e839-41fd-b878-aa0bcaf2fb16",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Min bottom elevation 44.10 m., max bottom elevation 67.55 m.\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([44.1 , 45.33, 46.56, 47.8 , 49.03, 50.27, 51.5 , 52.74, 53.97,\n",
       "       55.21, 56.44, 57.67, 58.91, 60.14, 61.38, 62.61, 63.85, 65.08,\n",
       "       66.32, 67.55])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# get raster minimum and maximum \n",
    "minElev = np.nanmin(lakeBottom)\n",
    "maxElev = np.nanmax(lakeBottom)\n",
    "print('Min bottom elevation %.2f m., max bottom elevation %.2f m.'%(minElev,maxElev))\n",
    "\n",
    "# steps for calculation\n",
    "nSteps = 20\n",
    "\n",
    "# lake bottom elevation intervals\n",
    "elevSteps = np.round(np.linspace(minElev,maxElev,nSteps),2)\n",
    "elevSteps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "db4a7ca4-e88e-44f3-af56-0702418169de",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# definition of volume function\n",
    "def calculateVol(elevStep,elevDem,lakeRst):\n",
    "    tempDem = elevStep - elevDem[elevDem<elevStep]\n",
    "    tempVol = tempDem.sum()*lakeRst.res[0]*lakeRst.res[1]\n",
    "    return tempVol"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "98572b92-3f8e-4e30-9644-2e91a23fa649",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lake bottom elevations [44.1  45.33 46.56 47.8  49.03 50.27 51.5  52.74 53.97 55.21 56.44 57.67\n",
      " 58.91 60.14 61.38 62.61 63.85 65.08 66.32 67.55]\n",
      "Lake volume in million of cubic meters [0.0, 0.0, 0.13, 0.38, 0.68, 1.02, 1.4, 1.83, 2.31, 2.88, 3.56, 4.3, 5.1, 5.97, 6.91, 7.92, 8.99, 10.12, 11.3, 12.54]\n"
     ]
    }
   ],
   "source": [
    "# calculate volumes for each elevation\n",
    "volArray = []\n",
    "for elev in elevSteps:\n",
    "    tempVol = calculateVol(elev,lakeBottom,lakeRst)\n",
    "    volArray.append(tempVol)\n",
    "    \n",
    "print(\"Lake bottom elevations %s\"%elevSteps)\n",
    "volArrayMCM = [round(i/1000000,2) for i in volArray]\n",
    "print(\"Lake volume in million of cubic meters %s\"%volArrayMCM)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "71952366-040f-4b36-b400-58c0149ac29d",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 20448 characters]\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot values\n",
    "fig, ax = plt.subplots(figsize=(12,5))\n",
    "ax.plot(volArrayMCM,elevSteps,label='Patillas lake')\n",
    "ax.grid()\n",
    "ax.legend()\n",
    "ax.set_xlabel('Volume MCM')\n",
    "ax.set_ylabel('Elevation (masl)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "58f201de-1811-42bf-adbd-435e10cd4359",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 143584 characters]\n",
      "text/plain": [
       "<Figure size 1440x576 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot values\n",
    "fig, [ax1, ax2] = plt.subplots(1,2,figsize=(20,8),gridspec_kw={'width_ratios': [2, 1]})\n",
    "ax1.set_title('Lake bottom elevation')\n",
    "botElev = ax1.imshow(lakeBottom)\n",
    "\n",
    "divider = make_axes_locatable(ax1)\n",
    "cax = divider.append_axes('bottom', size='5%', pad=0.5)\n",
    "fig.colorbar(botElev, cax=cax, orientation='horizontal', label='Elevation (masl)') \n",
    "\n",
    "ax2.plot(volArrayMCM,elevSteps,label='Patillas lake')\n",
    "ax2.grid()\n",
    "ax2.legend()\n",
    "ax2.set_xlabel('Volume MCM')\n",
    "ax2.set_ylabel('Elevation (masl)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5eda84d-2db4-4b4b-9508-5d6d18ccf28c",
   "metadata": {},
   "source": [
    "<hr/>\n",
    "\n",
    "### Comparison with the USGS survey\n",
    "\n",
    "From this pu"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "f3384479-df07-4c4b-a05c-25aa0ba96574",
   "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>OBJECTID</th>\n",
       "      <th>Pool_Elevation_m</th>\n",
       "      <th>Storage_capacity_mcm</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>1</td>\n",
       "      <td>67.55</td>\n",
       "      <td>12.96</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2</td>\n",
       "      <td>66.55</td>\n",
       "      <td>11.79</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>3</td>\n",
       "      <td>65.55</td>\n",
       "      <td>10.68</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>4</td>\n",
       "      <td>64.55</td>\n",
       "      <td>9.64</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>5</td>\n",
       "      <td>63.55</td>\n",
       "      <td>8.70</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "   OBJECTID  Pool_Elevation_m  Storage_capacity_mcm\n",
       "0         1             67.55                 12.96\n",
       "1         2             66.55                 11.79\n",
       "2         3             65.55                 10.68\n",
       "3         4             64.55                  9.64\n",
       "4         5             63.55                  8.70"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "patVol = pd.read_csv('../Txt/patillasVolume.csv')\n",
    "patVol.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "13ac774d-2caa-4ac4-8148-a1624e871ae0",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Elevation (masl)')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "[encoded data, 25088 characters]\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize=(12,5))\n",
    "ax.plot(volArrayMCM,elevSteps,label='Patillas lake',ls= '-.')\n",
    "ax.plot(patVol['Storage_capacity_mcm'],patVol['Pool_Elevation_m'],label='Patillas lake - USGS Survey')\n",
    "ax.grid()\n",
    "ax.legend()\n",
    "ax.set_xlabel('Volume MCM')\n",
    "ax.set_ylabel('Elevation (masl)')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
