python/how-tos/How_to_Extract_Points.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to: Extracting EMIT Spectra at Specified Coordinates\n",
    "\n",
    "**Summary**  \n",
    "\n",
    "In this notebook we will open a netCDF4 file from the Earth Surface Minteral Dust Source Investigation (EMIT) as an `xarray.Dataset`. We will then extract the spectra at point coordinates from a `.csv` as a dataframe, then save and plot the data.\n",
    "\n",
    "**Requirements:**\n",
    "+ A NASA [Earthdata Login](https://urs.earthdata.nasa.gov/) account is required to download EMIT data   \n",
    "+ Selected the `emit_tutorials` environment as the kernel for this notebook.\n",
    "  + For instructions on setting up the environment, follow the the `setup_instructions.md` included in the `/setup/` folder of the repository.  \n",
    "\n",
    "**Learning Objectives**\n",
    "+ How to open an EMIT file as an `xarray.Dataset`\n",
    "+ How to extract spectra at coordinates listed in a `.csv` file "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Import the required Python libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Packages\n",
    "import sys\n",
    "import os\n",
    "import earthaccess\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import hvplot.pandas\n",
    "import hvplot.xarray\n",
    "import holoviews as hv\n",
    "sys.path.append('../modules/')\n",
    "from emit_tools import emit_xarray"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Login to your NASA Earthdata account and create a `.netrc` file using the `login` function from the `earthaccess` library. If you do not have an Earthdata Account, you can create one [here](https://urs.earthdata.nasa.gov/home). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "earthaccess.login(persist=True)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this notebook we will download the files necessary using `earthaccess`. You can also access the data in place or stream it, but this can slow due to the file sizes. Provide a URL for an EMIT L2A Reflectance granule."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "url = 'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2ARFL.001/EMIT_L2A_RFL_001_20220903T163129_2224611_012/EMIT_L2A_RFL_001_20220903T163129_2224611_012.nc'"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get an HTTPS Session using your earthdata login, set a local path to save the file, and download the granule asset - This may take a while, the reflectance file is approximately 1.8 GB."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get requests https Session using Earthdata Login Info\n",
    "fs = earthaccess.get_requests_https_session()\n",
    "# Retrieve granule asset ID from URL (to maintain existing naming convention)\n",
    "granule_asset_id = url.split('/')[-1]\n",
    "# Define Local Filepath\n",
    "fp = f'../../data/{granule_asset_id}'\n",
    "# Download the Granule Asset if it doesn't exist\n",
    "if not os.path.isfile(fp):\n",
    "    with fs.get(url,stream=True) as src:\n",
    "        with open(fp,'wb') as dst:\n",
    "            for chunk in src.iter_content(chunk_size=64*1024*1024):\n",
    "                dst.write(chunk)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open the file downloaded and defined as `fp`. To do this, we will use the `emit_tools` module which contains a few helpful functions that can be used with EMIT data. Use the `ortho=True` option to orthorectify the dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = emit_xarray(fp, ortho=True)\n",
    "ds"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now open the `.csv` included in the `/data/` directory as a `pandas.dataframe`.\n",
    "\n",
    "> Note: The category values here are arbitrary and included as an example of an additional column users may want."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = pd.read_csv('../../data/sample_coords.csv')\n",
    "points"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make a plot to visualize the points we're going to select on the dataset. Here we use the reflectance values at 850nm as our basemap. To mask fill-values, assign them to `np.nan`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Assign fill-values to NaN\n",
    "ds.reflectance.data[ds.reflectance.data == -9999] = np.nan\n",
    "# Select a single band and plot points on top\n",
    "ds.sel(wavelengths=850,method='nearest').hvplot.image(cmap='greys', frame_height=600, frame_width=600, geo=True, crs='EPSG:4326').opts(title=f\"Example Points\")*\\\n",
    "points.hvplot.points(x='Longitude',y='Latitude', by='ID', color=hv.Cycle('Dark2'), geo=True, crs='EPSG:4326')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the `points` dataframe index as `ID` to utilize our existing point ID's as an index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = points.set_index(['ID'])"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Convert the dataframe to an `xarray.Dataset`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xp = points.to_xarray()\n",
    "xp"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Select the data from our EMIT dataset using the Latitude and Longitude coordinates from our point dataset, then convert the output to a `pandas.dataframe`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "extracted = ds.sel(latitude=xp.Latitude,longitude=xp.Longitude, method='nearest').to_dataframe()\n",
    "extracted"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output is a longform dataframe using the `'ID'` field as an index. This is missing our `'Category'` column from our original dataframe. Use the `pd.join` function to add the `'Category'` column to our dataset using `'ID'` as an index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = extracted.join(points['Category'], on=['ID'])\n",
    "df"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we have a dataframe containing our initial data, in addition to the extracted point data. This a a good place to save an output as a `.csv`. Go ahead and do that below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.to_csv('../../data/example_out.csv')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use our dataframe to plot the reflectance data we extracted, but first, mask the reflectance values of `-0.01`, which represent deep water vapor absorption regions. To do this, assign values where the reflectance = `-0.01` to `np.nan`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Data values are slightly different than -0.01 because of floating point precision, which requires a tolerance to be set\n",
    "df['reflectance'] = df['reflectance'].where(abs(df['reflectance'] - (-0.01)) > 1e-8, np.nan)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the data using `hvplot`. We can use `by=` to separate the reflectances by their `ID`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "df.hvplot(x='wavelengths',y='reflectance', by=['ID'], color=hv.Cycle('Dark2'), frame_height=400, frame_width=600).opts(title='Example Points - Reflectance', xlabel='Wavelengths (nm)',ylabel='Reflectance')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contact Info:  \n",
    "\n",
    "Email: LPDAAC@usgs.gov  \n",
    "Voice: +1-866-573-3222  \n",
    "Organization: Land Processes Distributed Active Archive Center (LP DAAC)¹  \n",
    "Website: <https://lpdaac.usgs.gov/>  \n",
    "Date last modified: 03-13-2024  \n",
    "\n",
    "¹Work performed under USGS contract G15PD00467 for NASA contract NNG14HH33I.  "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.13"
  },
  "vscode": {
   "interpreter": {
    "hash": "3292b2aceff7d39327a7519422d4180a7c9b133202090f26e797e3dd8f2c7877"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}