{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "22434761",
   "metadata": {},
   "source": [
    "# 3 Equivalent Water Thickness/Canopy Water Content from Imaging Spectroscopy Data\n",
    "\n",
    "> This notebook is part of a numbered series (01-05) that focuses specifically on combining EMIT VSWIR and ECOSTRESS TIR data for integrated land surface analysis.\n",
    "\n",
    "**Summary**  \n",
    "\n",
    "In this notebook we will explore how Equivalent Water Thickness (EWT) or Canopy Water Content (CWC) can be calculated from the Earth Surface Mineral Dust Source Investigation (EMIT) L2A Reflectance Product, then we will apply this knowledge to calculate CWC over the [Jack and Laura Dangermond Preserve](https://www.dangermondpreserve.org/) located near Santa Barbara, CA. \n",
    "\n",
    "<div>\n",
    "<img src=\"../../img/canopy_water_content_fig.png\" width=\"750\"/>\n",
    "</div>\n",
    "\n",
    "**Background**\n",
    "\n",
    "Equivalent Water Thickness (EWT) is the predicted thickness or absorption path length in centimeters (cm) of water that would be required to yield an observed spectra. In the context of vegetation, this is equivalent to canopy water content (CWC) in g/cm^2 because a cm^3 of water has a mass of 1g. \n",
    "\n",
    "CWC can be derived from surface reflectance spectra because they provide information about the composition of the target, including water content. Reflectance is the fraction of incoming solar radiation reflected by Earth's surface. Different materials reflect varying proportions of radiation based upon their chemical composition and physical properties, giving materials their own unique spectral signature or fingerprint. In particular, liquid water causes characteristic absorption features to appear in the near-infrared wavelengths of the solar spectrum, which enables an estimation of its content.\n",
    "\n",
    "CWC correlates with vegetation type and health, as well as wildfire risk. The methods used here to calculate CWC are based on the [ISOFIT python package](https://github.com/isofit/isofit/tree/main). The Beer-Lambert physical model used to calculate CWC is described in [Green et al. (2006)](https://doi.org/10.1029/2005WR004509) and [Bohn et al. (2020)](https://doi.org/10.1016/j.rse.2020.111708). It uses wavelength-dependent absorption coefficients of liquid water to determine the absorption path length as a function of absorption feature depth. Of note, this model does not account for multiple scattering effects within the canopy and may result in overestimation of CWC (Bohn et al., 2020).\n",
    "\n",
    "The [Jack and Laura Dangermond Preserve](https://www.dangermondpreserve.org/) and its surrounding lands are one of the last \"wild coastal\" regions in Southern California. The preserve is over 24,000 acres and consists of several ecosystem types and is home to over 600 plant species and over 200 wildlife species.\n",
    "\n",
    "More about the [EMIT mission](https://earth.jpl.nasa.gov/emit/) and [EMIT products](https://lpdaac.usgs.gov/product_search/?query=emit&status=Operational&view=cards&sort=title).\n",
    "\n",
    "**References**\n",
    "\n",
    "- Shrestha, Rupesh. 2023. Equivalent water thickness/canopy water content from hyperspectral data. Jupyter Notebook. Oak Ridge National Laboratory Distributed Active Archive Center. https://github.com/rupesh2/ewt_cwc/tree/main\n",
    "- Bohn, N., L. Guanter, T. Kuester, R. Preusker, and K. Segl. 2020. Coupled retrieval of the three phases of water from spaceborne imaging spectroscopy measurements. Remote Sensing of Environment 242:111708. https://doi.org/10.1016/j.rse.2020.111708\n",
    "- Green, R.O., T.H. Painter, D.A. Roberts, and J. Dozier. 2006. Measuring the expressed abundance of the three phases of water with an imaging spectrometer over melting snow. Water Resources Research 42:W10402. https://doi.org/10.1029/2005WR004509\n",
    "- Thompson, D.R., V. Natraj, R.O. Green, M.C. Helmlinger, B.-C. Gao, and M.L. Eastwood. 2018. Optimal estimation for imaging spectrometer atmospheric correction. Remote Sensing of Environment 216:355–373. https://doi.org/10.1016/j.rse.2018.07.003\n",
    "\n",
    "**Requirements** \n",
    " - [NASA Earthdata Account](https://urs.earthdata.nasa.gov/home)  \n",
    " - *No Python setup requirements if connected to the workshop cloud instance!*  \n",
    " - **Local Only** Set up Python Environment - See **setup_instructions.md** in the `/setup/` folder to set up a local compatible Python environment  \n",
    " - Downloaded necessary files. This is done at the end of the [01_Finding_Concurrent_Data](01_Finding_Concurrent_Data.ipynb) notebook.  \n",
    "\n",
    "**Learning Objectives**  \n",
    "- Calculate CWC of a single pixel  \n",
    "- Calculate CWC of an ROI  \n",
    "\n",
    "**Tutorial Outline**\n",
    "\n",
    "3.1 Setup  \n",
    "3.2 Opening EMIT Data  \n",
    "3.3 Extracting Reflectance of a Pixel  \n",
    "3.4 Calculating CWC   \n",
    "    3.4.1 Single Point  \n",
    "    3.4.2 DataFrame of Points  \n",
    "3.5 Applying Inversion in Parallel Across an ROI  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79781004",
   "metadata": {},
   "source": [
    "## 3.1 Setup\n",
    "\n",
    "Uncomment the below cell and install the `ray` python library. This is only necessary if using the 2i2c Cloud Instance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f2ecb342",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# !pip install \"ray[default]\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2d3b159",
   "metadata": {},
   "source": [
    "Import the required libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4b9dbb69",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Import Packages\n",
    "import os\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "from osgeo import gdal\n",
    "import rasterio as rio\n",
    "import rioxarray as rxr\n",
    "from matplotlib import pyplot as plt\n",
    "import hvplot.xarray\n",
    "import hvplot.pandas\n",
    "import pandas as pd\n",
    "import earthaccess\n",
    "\n",
    "from vitals.emit_tools import emit_xarray\n",
    "from vitals.ewt_calc import calc_ewt\n",
    "from scipy.optimize import least_squares"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "55d95b68",
   "metadata": {},
   "source": [
    "## 3.2 Open an EMIT File\n",
    "\n",
    "EMIT L2A Reflectance Data are distributed in a non-orthocorrected spatially raw NetCDF4 (.nc) format consisting of the data and its associated metadata. To work with this data, we will use the `emit_xarray` function from the `emit_tools.py` module included in the repository."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb74b0d7",
   "metadata": {},
   "source": [
    "### 3.2.1 Streaming Data\n",
    "\n",
    "As we did in the previous notebook, to stream the data, log in using `earthaccess`, then start and https session and open the url for the scene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "42a89f9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Login to NASA Earthdata\n",
    "earthaccess.login(persist=True)\n",
    "\n",
    "# Get Https Session using Earthdata Login Info\n",
    "fs = earthaccess.get_fsspec_https_session()\n",
    "\n",
    "# Define Local Filepath\n",
    "fp = fs.open('https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2ARFL.001/EMIT_L2A_RFL_001_20230401T203751_2309114_002/EMIT_L2A_RFL_001_20230401T203751_2309114_002.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a64718b",
   "metadata": {},
   "source": [
    "### 3.2.2 Downloading Data\n",
    "\n",
    "If you have downloaded the data, you can just set the filepath."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ac09f0a2",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# fp = '../../data/EMIT_L2A_RFL_001_20230401T203751_2309114_002.nc'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "da88f613",
   "metadata": {},
   "source": [
    "Open the file using the `emit_xarray` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a8abd337",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ds = emit_xarray(fp,ortho=True).load()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "302737ae",
   "metadata": {},
   "source": [
    "Set the fill values equal to `np.nan` improved visualization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8c866b27",
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.reflectance.data[ds.reflectance.data == -9999] = np.nan"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "849f0ea4",
   "metadata": {},
   "source": [
    "Similarly to what we did in the previous notebook, we can plot a single band to get an idea of where our scene is. This is the same scene we processed earlier."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "466fa65e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "emit_layer = ds.sel(wavelengths=850,method='nearest')\n",
    "emit_layer.hvplot.image(cmap='viridis',geo=True, tiles='ESRI', frame_width=720,frame_height=405, alpha=0.7, fontscale=2).opts(\n",
    "    title=f\"{emit_layer.wavelengths:.3f} {emit_layer.wavelengths.units}\", xlabel='Longitude',ylabel='Latitude')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "683bb1d0",
   "metadata": {},
   "source": [
    "## 3.3 Extracting Reflectance of a Single Pixel"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5756d7b9",
   "metadata": {},
   "source": [
    "We can mask out the -.01 values used to represent the region of the spectra with strong atmospheric water vapor absorption features."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "58667aca",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ds['reflectance'].data[:,:,ds['good_wavelengths'].data==0] = np.nan"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a759ad1e",
   "metadata": {},
   "source": [
    "Retrieve the spectra from a sample point by providing a latitude and longitude along with a method using the `sel` function. This will select the pixel closest to the provided coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f5d6fbbb",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "point = ds.sel(latitude=34.5399,longitude=-120.3529, method='nearest')\n",
    "point"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5aa525a5",
   "metadata": {},
   "source": [
    "We can plot this information to see the spectra."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a4c7327e",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "point.hvplot.line(x='wavelengths',y='reflectance',color='black').opts(title=f\"Latitude: {point.latitude.values:.3f} Longitude: {point.longitude.values:.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "21ffb183",
   "metadata": {},
   "source": [
    "## 3.4 Calculating CWC\n",
    "\n",
    "As mentioned in the background we use the surface reflectance to estimate CWC. The unique spectral signatures allow identification and quantification based upon the wavelength-dependent absorption coefficients of liquid water. The EMIT mission has applied similar approaches to identify dust source minerals as well as methane point source emissions. The path length of liquid water absorption can be estimated by utilizing a least squares inversion to minimize the residuals between the EMIT reflectance and the Beer-Lambert model (Green et al.,2006), which relates the wavelength-dependent absorption to the path length the photons are traveling through the material. During the inversion, the path lengths are iteratively adjusted to match the modeled spectra to the EMIT reflectance within the water absorption feature region from 850 to 1100 nm.\n",
    "\n",
    "First, define a Beer-Lambert Model function that returns the vector of residuals between measured and modeled surface reflectance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "78227ca3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# https://github.com/isofit/isofit/blob/main/isofit/inversion/inverse_simple.py#L514C1-L532C17\n",
    "def beer_lambert_model(x, y, wl, alpha_lw):\n",
    "    \"\"\"Function, which computes the vector of residuals between measured and modeled surface reflectance optimizing\n",
    "    for path length of surface liquid water based on the Beer-Lambert attenuation law.\n",
    "\n",
    "    Args:\n",
    "        x:        state vector (liquid water path length, intercept, slope)\n",
    "        y:        measurement (surface reflectance spectrum)\n",
    "        wl:       instrument wavelengths\n",
    "        alpha_lw: wavelength dependent absorption coefficients of liquid water\n",
    "\n",
    "    Returns:\n",
    "        resid: residual between modeled and measured surface reflectance\n",
    "    \"\"\"\n",
    "\n",
    "    attenuation = np.exp(-x[0] * 1e7 * alpha_lw)\n",
    "    rho = (x[1] + x[2] * wl) * attenuation\n",
    "    resid = rho - y\n",
    "\n",
    "    return resid"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe4a739e",
   "metadata": {},
   "source": [
    "We need some lab measurements of the complex refractive index of liquid water to obtain the wavelength-dependent absorption coefficients. They are calculated by taking four times the product of Pi and the imaginary part of the refractive index, divided by wavelength. The refractive index of liquid water per wavelength is provided by the `k_liquid_water_ice.csv` in the `data` folder. We can also preview this data to get a better understanding of the information we are using."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf91083c",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "wp_fp = '../../data/k_liquid_water_ice.csv'\n",
    "k_wi = pd.read_csv(wp_fp)\n",
    "k_wi.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0894db94",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(2,4, figsize=(15, 6),  sharex=True, sharey=True, constrained_layout=True)\n",
    "axs = axs.ravel()\n",
    "col_n = 0\n",
    "for i in range(0, 7):\n",
    "    x = k_wi.iloc[:, col_n+i]\n",
    "    y = k_wi.iloc[:, col_n+i+1]\n",
    "    axs[i].scatter(x, y)\n",
    "    axs[i].set_title(y.name)\n",
    "    col_n+=1\n",
    "fig.supylabel('imaginary parts of refractive index')\n",
    "fig.supxlabel('wavelength')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d89cd2c0",
   "metadata": {},
   "source": [
    "Now define a function to get the desired data from the csv file."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3928356b",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# https://github.com/isofit/isofit/blob/dev/isofit/core/common.py#L461C1-L488C26\n",
    "def get_refractive_index(k_wi, a, b, col_wvl, col_k):\n",
    "    \"\"\"Convert refractive index table entries to numpy array.\n",
    "\n",
    "    Args:\n",
    "        k_wi:    variable\n",
    "        a:       start line\n",
    "        b:       end line\n",
    "        col_wvl: wavelength column in pandas table\n",
    "        col_k:   k column in pandas table\n",
    "\n",
    "    Returns:\n",
    "        wvl_arr: array of wavelengths\n",
    "        k_arr:   array of imaginary parts of refractive index\n",
    "    \"\"\"\n",
    "\n",
    "    wvl_ = []\n",
    "    k_ = []\n",
    "\n",
    "    for ii in range(a, b):\n",
    "        wvl = k_wi.at[ii, col_wvl]\n",
    "        k = k_wi.at[ii, col_k]\n",
    "        wvl_.append(wvl)\n",
    "        k_.append(k)\n",
    "\n",
    "    wvl_arr = np.asarray(wvl_)\n",
    "    k_arr = np.asarray(k_)\n",
    "\n",
    "    return wvl_arr, k_arr"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fa3938a8",
   "metadata": {},
   "source": [
    "Lastly, to calculate CWC we define a function that uses least squares optimization to minimize the residuals of our Beer-Lambert Model and find a likely path length of liquid water."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9559b555",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# https://github.com/isofit/isofit/blob/main/isofit/inversion/inverse_simple.py#L443C1-L511C24\n",
    "def invert_liquid_water(\n",
    "    rfl_meas: np.array,\n",
    "    wl: np.array,\n",
    "    l_shoulder: float = 850,\n",
    "    r_shoulder: float = 1100,\n",
    "    lw_init: tuple = (0.02, 0.3, 0.0002),\n",
    "    lw_bounds: tuple = ([0, 0.5], [0, 1.0], [-0.0004, 0.0004]),\n",
    "    ewt_detection_limit: float = 0.5,\n",
    "    return_abs_co: bool = False,\n",
    "):\n",
    "    \"\"\"Given a reflectance estimate, fit a state vector including liquid water path length\n",
    "    based on a simple Beer-Lambert surface model.\n",
    "\n",
    "    Args:\n",
    "        rfl_meas:            surface reflectance spectrum\n",
    "        wl:                  instrument wavelengths, must be same size as rfl_meas\n",
    "        l_shoulder:          wavelength of left absorption feature shoulder\n",
    "        r_shoulder:          wavelength of right absorption feature shoulder\n",
    "        lw_init:             initial guess for liquid water path length, intercept, and slope\n",
    "        lw_bounds:           lower and upper bounds for liquid water path length, intercept, and slope\n",
    "        ewt_detection_limit: upper detection limit for ewt\n",
    "        return_abs_co:       if True, returns absorption coefficients of liquid water\n",
    "\n",
    "    Returns:\n",
    "        solution: estimated liquid water path length, intercept, and slope based on a given surface reflectance\n",
    "    \"\"\"\n",
    "    \n",
    "    # Ensure least squares is done with float64 datatype (added)\n",
    "    wl = np.float64(wl)\n",
    "    \n",
    "    # params needed for liquid water fitting\n",
    "    lw_feature_left = np.argmin(abs(l_shoulder - wl))\n",
    "    lw_feature_right = np.argmin(abs(r_shoulder - wl))\n",
    "    wl_sel = wl[lw_feature_left : lw_feature_right + 1]\n",
    "\n",
    "    # adjust upper detection limit for ewt if specified\n",
    "    if ewt_detection_limit != 0.5:\n",
    "        lw_bounds[0][1] = ewt_detection_limit\n",
    "\n",
    "    # load imaginary part of liquid water refractive index and calculate wavelength dependent absorption coefficient\n",
    "    # __file__ should live at isofit/isofit/inversion/\n",
    "    \n",
    "    \n",
    "    data_dir_path = \"../../data/\"\n",
    "    path_k = os.path.join(data_dir_path,\"k_liquid_water_ice.csv\")\n",
    "    \n",
    "    #isofit_path = os.path.dirname(os.path.dirname(os.path.dirname(__file__)))\n",
    "    #path_k = os.path.join(isofit_path, \"data\", \"iop\", \"k_liquid_water_ice.xlsx\")\n",
    "\n",
    "    # k_wi = pd.read_excel(io=path_k, sheet_name=\"Sheet1\", engine=\"openpyxl\")\n",
    "    # wl_water, k_water = get_refractive_index(\n",
    "    #     k_wi=k_wi, a=0, b=982, col_wvl=\"wvl_6\", col_k=\"T = 20°C\"\n",
    "    # )\n",
    "    k_wi = pd.read_csv(path_k)\n",
    "    wl_water, k_water = get_refractive_index(\n",
    "        k_wi=k_wi, a=0, b=982, col_wvl=\"wvl_6\", col_k=\"T = 20°C\"\n",
    "    )\n",
    "    kw = np.interp(x=wl_sel, xp=wl_water, fp=k_water)\n",
    "    abs_co_w = 4 * np.pi * kw / wl_sel\n",
    "\n",
    "    rfl_meas_sel = rfl_meas[lw_feature_left : lw_feature_right + 1]\n",
    "\n",
    "    x_opt = least_squares(\n",
    "        fun=beer_lambert_model,\n",
    "        x0=lw_init,\n",
    "        jac=\"2-point\",\n",
    "        method=\"trf\",\n",
    "        bounds=(\n",
    "            np.array([lw_bounds[ii][0] for ii in range(3)]),\n",
    "            np.array([lw_bounds[ii][1] for ii in range(3)]),\n",
    "        ),\n",
    "        max_nfev=15,\n",
    "        args=(rfl_meas_sel, wl_sel, abs_co_w),\n",
    "    )\n",
    "\n",
    "    solution = x_opt.x\n",
    "\n",
    "    if return_abs_co:\n",
    "        return solution, abs_co_w\n",
    "    else:\n",
    "        return solution"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f55c0db9",
   "metadata": {},
   "source": [
    "### 3.4.1 CWC of a Single Point\n",
    "\n",
    "Now that we have all of the pieces, we can estimate the CWC of a single pixel using our `invert_liquid_water` function. By default there is a detection limit of 0.5, if we are hitting this threshold we can adjust by changing the `ewt_detection_limit` argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67e03abe",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ewt = invert_liquid_water(point.reflectance.values,point.wavelengths.values)\n",
    "print(f\"EWT for ({point.longitude.values:.3f},{point.latitude.values:.3f}): {ewt[0]:.3f} cm\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2494446f",
   "metadata": {},
   "source": [
    "### 3.4.2 CWC of a DataFrame of Points\n",
    "\n",
    "We can also apply this function to the point data we selected in the previous notebook with the interactive plot.\n",
    "\n",
    "Read in the csv file we created."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de554c8f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "points_df = pd.read_csv(\"../../data/emit_dangermond_click_data.csv\")\n",
    "points_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c155fea",
   "metadata": {},
   "source": [
    "Now create an array of wavelengths from the column names starting with the first wavelength value (column 3)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ff46c91d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Get wavelength values\n",
    "wavelength_values = points_df.columns[3::].to_numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e6b8ad8",
   "metadata": {
    "tags": []
   },
   "source": [
    "Iterate by row through our dataframe, selecting the reflectance values and providing them to the `invert_liquid_water` function. Afterwards, add an CWC column to our dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0d111b14",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "# Create empty list\n",
    "ewt_values = []\n",
    "# Iterate through rows\n",
    "for _i in points_df.index.to_list():\n",
    "    # Get reflectance array to pass to function\n",
    "    rfl_values = points_df.iloc[_i,3::]\n",
    "    # Use invert liquid water function and append results to list\n",
    "    ewt_values.append(invert_liquid_water(rfl_values,wavelength_values)[0])    \n",
    "# Add to our existing dataframe at Column Index 3\n",
    "points_df.insert(3, \"ewt\", ewt_values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60ca9bfc",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "points_df"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "24e3f70c",
   "metadata": {},
   "source": [
    "For larger sets of point data we would want to do this in parallel. We can also calculate CWC for an area, where we definitely need to utilize parallel processing."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94af67c3",
   "metadata": {},
   "source": [
    "## 3.5 CWC Calculation of an ROI\n",
    "\n",
    "In the previous notebook, we subset our region of interest and exported the file. Since the CWC calculation is computationally intensive, it can take a while to process large scenes, so it is more efficient to do this spatial subsetting up front. We can use a function included in the `ewt_calc.py` module to calculate CWC on a cropped image, and create a cloud-optimized GeoTIFF (COG) file containing the results.\n",
    "\n",
    "Set our input filepaths and output directory. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "60225628",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "fp = '../../data/EMIT_L2A_RFL_001_20230401T203751_2309114_002_dangermond.nc'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4639996d",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "out_dir = '../../data/'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "34ee08dd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "roi_ds = xr.open_dataset(fp, decode_coords='all')\n",
    "roi_ds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b77d9888",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "roi_ds.sel(wavelengths=850,method='nearest').hvplot.image(x='longitude',y='latitude',cmap='viridis',geo=True, tiles='ESRI', frame_width=720,frame_height=405, alpha=0.7, fontscale=2).opts(\n",
    "    title=f\"Dangermond ROI - RFL at 850 nm\", xlabel='Longitude',ylabel='Latitude')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c597b8d3",
   "metadata": {},
   "source": [
    "Use the `calc_ewt` function to calculate CWC of the cropped image. This function will also create a COG file containing the CWC results. We can also specify the number of CPUs to use manually with a `n_cpu` argument, or leave it blank to use all but one of the available CPUs. If we set the `return_cwc` argument to true, the function will also return the COG. \n",
    "\n",
    "This will take some time, about 5 minutes, because we're doing the calculation for roughly 63,000 pixels. Also note that here we provide the `ewt_detection_limit` to increase it from the default of 0.5 in the function. We do this because there are several regions containing plants that hold significant quantities of water in this scene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b8cb40a8",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "ds_cwc = calc_ewt(fp, out_dir, ewt_detection_limit=1.5, return_cwc=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6a3b670f",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ds_cwc"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0f46c5f7",
   "metadata": {},
   "source": [
    "Plot CWC of our ROI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "03af8902",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "ds_cwc.hvplot.image(x='longitude',y='latitude',cmap='viridis',geo=True, tiles='ESRI', frame_width=720,frame_height=405, alpha=0.7, fontscale=2).opts(\n",
    "    title=f\"{ds_cwc.cwc.long_name} ({ds_cwc.cwc.units})\", xlabel='Longitude',ylabel='Latitude')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ea2964a",
   "metadata": {},
   "source": [
    "## Contact Info:\n",
    "Email: LPDAAC@usgs.gov  \n",
    "Voice: +1-866-573-3222  \n",
    "Organization: Land Processes Distributed Active Archive Center (LP DAAC)¹  \n",
    "Website: <https://www.earthdata.nasa.gov/centers/lp-daac>\n",
    "\n",
    "¹Work performed under USGS contract 140G0126D0001 for NASA contract NNG14HH33I."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "lp",
   "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.12.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}