python/how-tos/How_to_use_EMIT_Quality_data.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to: Use EMIT Quality Data\n",
    "\n",
    "There are two quality layers associated with L2A Reflectance Product contained within the L2A Mask file. The `mask` variable contains 6 binary flag bands (1-5, 8) which should be excluded from analysis, and two data bands (6,7).\n",
    "\n",
    "The second `band_mask` variable indicates whether or not any given wavelength of any given pixel is interpolated. Interpolation occurs either due to a focal plane array bad pixel, or from saturation. This data is provided as a packed unsigned integer array with 36 elements.\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 build a mask using the quality flags from an EMIT L2A Mask file\n",
    "+ How to build a mask by unpacking the `band_mask` data\n",
    "+ How to apply the band and quality masks to another EMIT file."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Using EMIT Quality Flag Data\n",
    "Import the required Python libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import Packages\n",
    "import os\n",
    "import earthaccess\n",
    "import netCDF4 as nc\n",
    "from osgeo import gdal\n",
    "import numpy as np\n",
    "import xarray as xr\n",
    "import hvplot.xarray\n",
    "import holoviews as hv\n",
    "import sys\n",
    "sys.path.append('../modules/')\n",
    "import emit_tools"
   ]
  },
  {
   "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)"
   ]
  },
  {
   "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 URLs for an EMIT L2A Reflectance and L2A Mask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "urls = ['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',\n",
    "        'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2ARFL.001/EMIT_L2A_RFL_001_20220903T163129_2224611_012/EMIT_L2A_MASK_001_20220903T163129_2224611_012.nc']\n",
    "# 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",
    "for url in urls:\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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set a filepath for the reflectance and mask files downloaded."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fp = '../../data/EMIT_L2A_RFL_001_20220903T163129_2224611_012.nc'\n",
    "fp_mask = '../../data/EMIT_L2A_MASK_001_20220903T163129_2224611_012.nc'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The most efficient way to utilize the mask is to apply it before orthorectification because the orthorectified datasets take up more space. To apply a  mask using the L2A Mask file, we want to open it, specify which bands to use in construction of a mask, and then apply the mask.\n",
    "\n",
    "To do this, first take a look at what each band will mask by reading in the `sensor_band_parameters` group from the mask file as an `xarray.dataset` then converting to a dataframe.\n",
    "\n",
    ">Note: In the user guide, the bands are indexed as 1-8 not 0-7 as used here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_parameters_ds = xr.open_dataset(fp_mask,engine = 'h5netcdf', group='sensor_band_parameters')\n",
    "mask_key = mask_parameters_ds['mask_bands'].to_dataframe()\n",
    "mask_key"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Visualizing the mask"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's start by reviewing some of the content of these bands, and seeing what's available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import matplotlib.gridspec as gridspec\n",
    "fig = plt.figure(figsize=(20,50))\n",
    "gs = gridspec.GridSpec(ncols=3, nrows=len(mask_key), figure=fig)\n",
    "\n",
    "ds = emit_tools.emit_xarray(fp, ortho = False)\n",
    "mask_ds = emit_tools.emit_xarray(fp_mask, ortho=False)\n",
    "\n",
    "rgb_inds = np.array([np.nanargmin(abs(ds['wavelengths'].values - x)) for x in [650, 560, 470]])\n",
    "rgb = ds['reflectance'].values[:,:,rgb_inds] # subset RGB\n",
    "rgb[rgb < 0] = np.nan\n",
    "rgb -= np.nanpercentile(rgb,2,axis=(0,1))[np.newaxis,np.newaxis,:] # scale from 2-95 %\n",
    "rgb /= np.nanpercentile(rgb,95,axis=(0,1))[np.newaxis,np.newaxis,:]\n",
    "\n",
    "for _n in range(int(len(mask_key)/2)):\n",
    "    ax = fig.add_subplot(gs[_n, 0])\n",
    "    plt.imshow(rgb);\n",
    "    plt.axis('off')\n",
    "    plt.title('RGB')\n",
    "    \n",
    "    ax = fig.add_subplot(gs[_n, 1])\n",
    "    md = mask_ds['mask'].values[...,_n]\n",
    "    md[np.isnan(rgb[...,0])] = np.nan\n",
    "    plt.imshow(md);\n",
    "    plt.axis('off')\n",
    "    plt.title(mask_key['mask_bands'][_n])\n",
    "    \n",
    "    ax = fig.add_subplot(gs[_n, 2])\n",
    "    md = mask_ds['mask'].values[...,_n+int(len(mask_key)/2)]\n",
    "    md[np.isnan(rgb[...,0])] = np.nan\n",
    "    plt.imshow(md);\n",
    "    plt.axis('off')\n",
    "    plt.title(mask_key['mask_bands'][_n+int(len(mask_key)/2)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some of these bands are direct masks (Cloud, Dilated, Currus, Water, Spacecraft), and some (AOD550 and H2O (g cm-2)) are information calculated during the L2A reflectance retrieval that may be used as additional screening, depending on the application.  The final mask that the EMIT mission will use for its minerological applications is shown as the Aggreged Flag - but not all users might want this particular mask."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a closer look at some of those bands with additional information, that could be used either to screen specific content out or that might be used as signal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "md = mask_ds['mask'].values[:,:,list(mask_key['mask_bands']).index('AOD550')]\n",
    "md[np.isnan(rgb[...,0])] = np.nan\n",
    "plt.imshow(md, vmin=np.nanpercentile(md,2),vmax=np.nanpercentile(md,98));\n",
    "plt.title('AOD550')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "md = mask_ds['mask'].values[:,:,list(mask_key['mask_bands']).index('H2O (g cm-2)')]\n",
    "md[np.isnan(rgb[...,0])] = np.nan\n",
    "plt.imshow(md, vmin=np.nanpercentile(md,2),vmax=np.nanpercentile(md,98));\n",
    "plt.title('H2O (g cm-2)')\n",
    "plt.colorbar()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Applying the Mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask_key"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The above dataframe shows exactly what each band contained within the file represents/will mask. For this example we will use flags 0,1,3, and 4 to remove any potential clouds and any artefacts caused by the space station. This can be done using the `quality_mask` function from the `emit_tools` module. This function combines the requested flags into a single mask and returns it as an array.\n",
    "\n",
    "Select the bands to use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flags = [0,1,3]\n",
    "flags"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now create the mask using the `quality_mask` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mask = emit_tools.quality_mask(fp_mask,flags)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see the regions of the unorthocorrected image that will be masked we can plot the mask array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flags = [0,1,3,4]\n",
    "mask = emit_tools.quality_mask(fp_mask,flags)\n",
    "fig = plt.figure(figsize=(15,15))\n",
    "gs = gridspec.GridSpec(ncols=2, nrows=2, figure=fig)\n",
    "\n",
    "ax = fig.add_subplot(gs[0, 0])\n",
    "plt.imshow(rgb)\n",
    "plt.scatter(1200,1200,c='red',marker='+')\n",
    "\n",
    "ax = fig.add_subplot(gs[0, 1])\n",
    "plt.imshow(mask)\n",
    "plt.scatter(1200,1200,c='red',marker='+')\n",
    "\n",
    "ax = fig.add_subplot(gs[1, :])\n",
    "plt.plot(ds['wavelengths'],ds['reflectance'].values[1200,1200,:])\n",
    "plt.xlabel('Wavelengths [nm]')\n",
    "plt.ylabel('Reflectance')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have a mask to apply, we can use it as the `qmask` parameter in the `emit_xarray` function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = emit_tools.emit_xarray(fp, ortho=True, qmask=mask)\n",
    "ds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This `ds` dataset is now orthorectified and the specified mask has been applied. Visualize the output using a plot of a red band (650 nm)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.reflectance.data[ds.reflectance.data == -9999] = np.nan\n",
    "ds.sel(wavelengths=650, method='nearest').hvplot.image(cmap='viridis', aspect = 'equal', frame_width=500, rasterize=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The modified `xarray.Dataset` can also be saved as a netCDF4 output that can be reopened using the `xarray.open_dataset` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds.to_netcdf('../../data/example_quality_nc_out.nc')\n",
    "# Example for Opening \n",
    "# ds = xr.open_dataset('../data/example_quality_nc_out.nc')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Using EMIT Band Mask Data\n",
    "\n",
    "The EMIT L2A Mask file also contains `band_mask` data, which indicates whether or not any given wavelength of any given pixel is interpolated. Interpolation occurs either due to a focal plane array bad pixel, or from saturation. This data comes as a packed unsigned integer array with 36 elements.\n",
    "\n",
    "Unpack the data an using the `band_mask` function from the `emit_tools` module. This function will unpack the data and create an array that can be used to mask the bands/pixels when added as an input into the `emit_xarray` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bmask = emit_tools.band_mask(fp_mask)\n",
    "bmask.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can quickly plot an example of a band where some pixels have been interpolated (band 234)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "plt.imshow(bmask[:,:,234])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or more helpfully, plot a representation of the crosstrack of the detector array to see where interpolation is occuring. If your research depends on spectral features contained within these interpolated crosstrack region, you may want to mask them out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(20,10))\n",
    "\n",
    "plt.imshow(bmask[0,...].T)\n",
    "plt.xlabel('Crosstrack Element')"
   ]
  },
  {
   "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: 07-03-2023  \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
}