python/tutorials/Generating_Methane_Spectral_Fingerprint.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generating Methane Spectral Fingerprint\n",
    "\n",
    "**Summary**  \n",
    "\n",
    "In this notebook, we'll examine a EMIT L1B At-Sensor Calibrated Radiance ([EMITL1BRAD](https://lpdaac.usgs.gov/products/emitl1bradv001/)) scene to visualize the spectra for a sample point inside a methane plume, and two sample points outside of methane plumes. To highlight the spectral differences we will look at an in-plume divided by out-of-plume ratio to get an example of the spectral signature or fingerprint of methane and visually compare this to the target methane signature.\n",
    "\n",
    "**Background**\n",
    "\n",
    "Methane is a major contributor to atmospheric radiative forcing because its more efficient at trapping radiation than other greenhouse gases, like carbon dioxide. Because the lifetime of methane in the atmosphere is only about 10 years, reducing methane emissions offers an effective way to curb anthropogenic contributions to atmospheric radiative forcing.\n",
    "\n",
    "The EMIT instrument is an imaging spectrometer that measures light in visible and infrared wavelengths. These measurements display unique spectral signatures that correspond to chemical composition. Although the primary mission focuses on mapping the mineral composition of Earth's surface, the EMIT instrument has also been used to successfully map methane point source emissions within its target mask using these unique spectral signatures. More details about EMIT and its associated products can be found in the **README.md** and on the [EMIT website](https://earth.jpl.nasa.gov/emit/).\n",
    "\n",
    "The EMITL1BRAD product provides at-sensor calibrated radiance values along with observation data in a spatially raw, non-orthocorrected format. This product is used to estimate the EMIT L2B Methane Enhancement Data ([EMITL2BCH4ENH](https://lpdaac.usgs.gov/products/emitl2bch4enhv001/)) product. The EMITL2BCH4ENH product is created using an adaptive matched filter to search each pixel's radiance spectrum for deviations that are characteristic of methane's absorption spectrum, and is only delivered to the LP DAAC for scenes where plume complexes have been identified. To reduce the risk of false positives, all EMITL2BCH4ENH data undergo a manual review process before being designated as a plume complex. For more information on the identification and manual review process, see Section 4.2.2 of the [EMIT GHG Algorithm Theoretical Basis Document (ATBD)](https://lpdaac.usgs.gov/documents/1696/EMIT_GHG_ATBD_V1.pdf). \n",
    "\n",
    "\n",
    "**References**\n",
    "\n",
    "Andrew K. Thorpe et al., Attribution of individual methane and carbon dioxide emission sources using EMIT observations from space. Sci. Adv.9, eadh2391 (2023). DOI:[10.1126/sciadv.adh2391](https://www.science.org/doi/10.1126/sciadv.adh2391)\n",
    "\n",
    "**Requirements** \n",
    " - Set up Python Environment - See **setup_instructions.md** in the `/setup/` folder\n",
    " - NASA Earthdata Login Account. [Sign Up](https://urs.earthdata.nasa.gov/users/new)\n",
    "\n",
    "**Data Used** \n",
    " - [EMIT L1B At-Sensor Calibrated Radiance and Geolocation Data (EMITL1BRAD)](https://lpdaac.usgs.gov/products/emitl1bradv001/)\n",
    " - [EMIT L2B Methane Enhancement Data (EMITL2BCH4ENH)](https://lpdaac.usgs.gov/products/emitl2arflv001/)\n",
    "\n",
    "\n",
    "**Learning Objectives** \n",
    " - Open and orthorectify an EMITL1BRAD scene\n",
    " - Open and visualize the scene and point data\n",
    " - Extract point data using coordinates in a .csv file\n",
    " - Visualize spectral fingerprint of methane and compare it to the modeled methane signature\n",
    "\n",
    "**Tutorial Outline**  \n",
    "\n",
    "1. [**Setup**](#setup)\n",
    "2. [**Opening EMIT Data**](#opening)\n",
    "3. [**Extracting Point Data**](#extracting_points)\n",
    "4. [**Methane Spectral Signature**](#spectral_signature)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Setup<a id='setup'></a>\n",
    "\n",
    "Import the necessary Python libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from osgeo import gdal\n",
    "import earthaccess\n",
    "import rasterio as rio\n",
    "import rioxarray as rxr\n",
    "import holoviews as hv\n",
    "import hvplot\n",
    "import hvplot.xarray\n",
    "import hvplot.pandas\n",
    "from skimage import exposure\n",
    "\n",
    "sys.path.append('../modules/')\n",
    "from emit_tools import emit_xarray"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download the EMIT L1B Radiance and L2B Methane Enhancement products for the scene we're going to look at."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "urls = ['https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL1BRAD.001/EMIT_L1B_RAD_001_20220815T042838_2222703_003/EMIT_L1B_RAD_001_20220815T042838_2222703_003.nc',\n",
    "        'https://data.lpdaac.earthdatacloud.nasa.gov/lp-prod-protected/EMITL2BCH4ENH.001/EMIT_L2B_CH4ENH_001_20220815T042838_2222703_003/EMIT_L2B_CH4ENH_001_20220815T042838_2222703_003.tif']\n",
    "\n",
    "# Authenticate and create an https session\n",
    "earthaccess.login(persist=True)\n",
    "fs = earthaccess.get_requests_https_session()\n",
    "\n",
    "for url in urls:\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)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the filepath for the radiance and methane enhancement files."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rad_fp = '../../data/EMIT_L1B_RAD_001_20220815T042838_2222703_003.nc'\n",
    "enh_fp = '../../data/EMIT_L2B_CH4ENH_001_20220815T042838_2222703_003.tif'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Opening EMIT Data<a id='opening'></a>\n",
    "\n",
    "The EMIT L1B At-Sensor Radiance data is distributed in a non-orthorectified spatially raw netCDF4 (.nc) format consisting of the data and its associated metadata. Inside the L1B file, there are 3 groups. \n",
    "\n",
    "1. The root group that can be considered the main dataset contains the radiance data described by the downtrack, crosstrack, and bands dimensions.  \n",
    "2. The sensor_band_parameters group containing the wavelength center and the full-width half maximum (FWHM) of each band.  \n",
    "3. The location group contains latitude and longitude values at the center of each pixel described by the crosstrack and downtrack dimensions, as well as a geometry lookup table (GLT) described by the ortho_x and ortho_y dimensions. The GLT is an orthorectified image (EPSG:4326) consisting of 2 layers containing downtrack and crosstrack indices. These index positions allow us to quickly project the raw data onto this geographic grid.  \n",
    "\n",
    "This data can be opened using the `netCDF4` and `xarray` libraries, or utilizing functions within the `emit_tools.py` module to organize them into a flattened (no groups) `xarray.Dataset` object. For this notebook, we will use functions from `emit_tools.py` to simplify working with the data. For more about the structure and use of `netCDF4` and `xarray` please see the [Exploring EMIT L2A Reflectance Jupyter Notebook](https://github.com/nasa/EMIT-Data-Resources/blob/main/python/tutorials/Exploring_EMIT_L2A_Reflectance.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open the radiance file using the `emit_xarray` function from the `emit_tools.py` module and orthorecitify it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rad = emit_xarray(rad_fp, ortho=True)\n",
    "rad"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The EMIT L2B Methane Enhancement Data represent an enhancement above background methane concentration for a 1 meter layer in parts-per-million (ppm) meter (m). These units are used rather than ppm because we are unable to measure the vertical extent of plumes. This data is distributed as a single band cloud-optimized geotiff (COG) and it has been orthocorrected. We can open this using the `rioxarray` library to place the data in an `xarray.DataArray` object."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Open the methane enhancement geotiff file using `rioxarray` and squeeze the `band` dimension to remove the extra dimension so our array will only have 2 dimensions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "enh = rxr.open_rasterio(enh_fp).squeeze('band',drop=True) \n",
    "enh"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Extracting Point Data <a id='extracting_points'></a>\n",
    "\n",
    "Open the `.csv` file included in the `/data/` directory as a `pandas.DataFrame`. This file contains the latitude and longitude coordinates for three points of interest. Two that are outside of a methane plume, and one that is inside."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define our Points for In-plume and out-of-plume\n",
    "points = pd.read_csv('../../data/methane_tutorial/methane_inout_points.csv')\n",
    "points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set the index as the `ID` column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "points = points.set_index(['ID'])\n",
    "points"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use the `sel` function from `xarray` to extract the radiance data at each point in our dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Extract target spectra from dataset\n",
    "point_ds = rad.sel(latitude=points.to_xarray().latitude, longitude=points.to_xarray().longitude, method='nearest')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After extracting, we can convert the data to a `pandas.DataFrame` and join it with our 'in-plume' column from the original dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "point_df = point_ds.to_dataframe().join(points['in-plume'],on=['ID'])\n",
    "point_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "At this point, we can save the data to a `.csv` file for future use."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# point_df.to_csv('../../data/methane_tutorial/point_df.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, visualize these points on an rgb image generated from the radiance file, and the methane enhancement data, just to get a better idea of where these points are located.\n",
    "\n",
    "To do this, first create an RGB data array from the radiance file using the `sel` function to select the bands nearest to the desired wavelengths."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an RGB from Radiance\n",
    "rgb = rad.sel(wavelengths=[650,560,470], method='nearest')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, use a function to rescale the brightness, so this image is easier to see."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "rgb.radiance.data[rgb.radiance.data == -9999] = 0\n",
    "rgb.radiance.data = exposure.rescale_intensity(rgb.radiance.data, in_range='image', out_range=(0,1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the necessary pieces to visualize we can make some spatial visualizations using `hvplot`. For the enhancement data, set the fill value of -9999 to `np.nan` to make it transparent."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create RGB Plot\n",
    "rgb_map = rgb.hvplot.rgb(x='longitude',y='latitude',bands='wavelengths',title='RGB Radiance', geo=True, crs='EPSG:4326')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create Methane Enhancement Plot\n",
    "enh.data[enh.data == -9999] = np.nan\n",
    "methane_map = enh.hvplot.image(x='x',y='y',cmap='viridis', geo=True, crs='EPSG:4326', clim=(0,enh.data.max()), title='Methane Enhancement', clabel='ppm m', xlabel='Longitude', ylabel='Latitude')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "point_map = point_df.hvplot.points(x='longitude',y='latitude', color='in-plume', cmap='HighContrast', geo=True, crs='EPSG:4326', hover=False, colorbar=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can combine these plots with an `*` operator to overlay them on the same plot. This does require that all are in the same CRS to display properly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(rgb_map*point_map) + (methane_map*point_map)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now lets look at the spectra. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "point_df.hvplot.line(x='wavelengths',y='radiance', by=['ID'], color=hv.Cycle('Dark2'), frame_height=400, frame_width=600, title = 'Radiance Spectra, ID 0 is in-plume' , xlabel='Wavelength (nm)', ylabel='Radiance (W/m^2/sr/nm)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Methane Spectral Signature<a id='spectral_signature'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's open a file containing the modeled methane signature and visualize it. This is the spectral fingerprint that we are looking for in EMIT data to identify methane plumes. Open our absorption coefficient file using `pandas`, add some column names and set an index."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open file \n",
    "ch4_ac = pd.read_csv('../../data/methane_tutorial/emit20220815t042838_ch4_target', sep='\\s+', header=None)\n",
    "# Add Column Names\n",
    "ch4_ac.columns = ['index','wavelength','value']\n",
    "# Set Index\n",
    "ch4_ac.set_index('index', inplace=True)\n",
    "ch4_ac"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a figure for the absorption coefficient. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ac_plot = ch4_ac.hvplot(x='wavelength',y='value', frame_height=400, frame_width=400, line_color='black', line_width=2, xlim=(2150,2450), ylim=(-1.5,0), xlabel='Wavelength (nm)', title='Methane Absorption Coefficient', ylabel='')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ac_plot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visualize a similar curve by creating a band ratio of in-plume to out-of-plume spectra. dividing the radiance for a pixel with methane by radiance for a pixel outside the plume to generate a diagnostic spectral fingerprint. This spectral fingerprint is confirmation that EMIT is observing methane enhancements with characteristic methane absorption features, which agree well with the modeled methane signature.\n",
    "\n",
    "This agreement is strong in the example we selected due to the similarity of the spectra without the contribution of methane. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create the band ratio, first separate our data into in-plume and out-of-plume dataframes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "in_plume = point_df.loc[point_df['in-plume'] == 1].copy()\n",
    "out_plume = point_df.loc[point_df['in-plume'] == 0].copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, add a column for the in/out band ratio using the in-plume divided by out-of-plume radiance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_plume['band_ratio'] = (in_plume.loc[0,'radiance']/out_plume['radiance'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create an `hvplot` object for our band ratio."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "in_out_plot = out_plume.hvplot(x='wavelengths',y='band_ratio', by=['ID'], color=hv.Cycle('Dark2'), frame_height=400, frame_width=400, xlim=(2150,2450), ylim=(0.85,1.05), ylabel='In Plume/Out of Plume Ratio', xlabel='Wavelength (nm)', title='In Plume/Out of Plume Ratio')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Overlay our absorption coefficient to show similarity between the two within the 2150 and 2450 nanometers (nm) range where methane spectral features are present. We can do this by setting our figure `xlim`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bokeh.models import GlyphRenderer, LinearAxis, LinearScale, Range1d\n",
    "\n",
    "def overlay_hook(plot, element):\n",
    "    # Adds right y-axis\n",
    "    p = plot.handles[\"plot\"]\n",
    "    p.extra_y_scales = {\"right\": LinearScale()}\n",
    "    p.extra_y_ranges = {\"right\": Range1d(-1.5,0)}\n",
    "    p.add_layout(LinearAxis(y_range_name=\"right\"), \"right\")\n",
    "\n",
    "   # find the last line and set it to right\n",
    "    lines = [p for p in p.renderers if isinstance(p, GlyphRenderer)]\n",
    "    lines[-1].y_range_name = \"right\"\n",
    "\n",
    "# Create Figure\n",
    "(in_out_plot.opts(ylim=(0.85,0.95)) * ac_plot.opts(color=\"k\")).opts(hooks=[overlay_hook]).opts(title='In Plume/Out of Plume and Absorption Coefficient') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can contrast this with our out-of-plume/out-of-plume band ratios to show how similar the in-plume and out-of-plume ratio is to this spectral signature. Create a column in our dataframe for the out-of-plume/out-of-plume band ratio and visualize it overlayed with our methane absorption feature using `hvplot`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_plume['out-out'] = (out_plume['radiance'] / out_plume.loc[2,'radiance'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "out_out_plot = out_plume.hvplot(x='wavelengths',y='out-out', by=['ID'], color=hv.Cycle('Dark2'), frame_height=400, frame_width=400, xlim=(2150,2450), ylim=(0.85,1.05), ylabel='Out/Out 2 Ratio', xlabel='Wavelength (nm)', title='Out of Plume/Out of Plume 2 Ratio')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from bokeh.models import GlyphRenderer, LinearAxis, LinearScale, Range1d\n",
    "\n",
    "def overlay_hook(plot, element):\n",
    "    # Adds right y-axis\n",
    "    p = plot.handles[\"plot\"]\n",
    "    p.extra_y_scales = {\"right\": LinearScale()}\n",
    "    p.extra_y_ranges = {\"right\": Range1d(-1.5,0)}\n",
    "    p.add_layout(LinearAxis(y_range_name=\"right\"), \"right\")\n",
    "\n",
    "   # find the last line and set it to right\n",
    "    lines = [p for p in p.renderers if isinstance(p, GlyphRenderer)]\n",
    "    lines[-1].y_range_name = \"right\"\n",
    "\n",
    "# Create Figure\n",
    "(out_out_plot * ac_plot.opts(color=\"k\")).opts(hooks=[overlay_hook]).opts(title='Out of Plume/Out of Plume and Absorption Coefficient') "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This highlights the similarity of the in-plume/out-of-plume band ratio to the modeled methane signature, as opposed to a band ratio where methane is not present. "
   ]
  },
  {
   "cell_type": "markdown",
   "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://lpdaac.usgs.gov/  \n",
    "Date last modified: 03-13-2024  \n",
    "\n",
    "¹Work performed under USGS contract G15PD00467 for NASA contract NNG14HH33I.  \n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "lpdaac_vitals_test",
   "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}