python/tutorials/GEDI_L1B_V2_Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Getting Started with GEDI L1B Version 2 Data in Python\n",
    "### This tutorial demonstrates how to work with the Geolocated Waveform ([GEDI01_B.002](https://doi.org/10.5067/GEDI/GEDI01_B.002)) data product.\n",
    "The Global Ecosystem Dynamics Investigation ([GEDI](https://lpdaac.usgs.gov/data/get-started-data/collection-overview/missions/gedi-overview/)) mission aims to characterize ecosystem structure and dynamics to enable radically improved quantification and understanding of the Earth's carbon cycle and biodiversity. The GEDI instrument produces high resolution laser ranging observations of the 3-dimensional structure of the Earth. GEDI is attached to the International Space Station and collects data globally between 51.6 N and 51.6 S latitudes at the highest resolution and densest sampling of any light detection and ranging (lidar) instrument in orbit to date. The Land Processes Distributed Active Archive Center (LP DAAC) distributes the GEDI Level 1 and Level 2 Version 1 and Version 2 products. The L1B and L2 GEDI products are archived and distributed in the HDF-EOS5 file format. \n",
    "\n",
    "---\n",
    "## Use Case Example:  \n",
    "This tutorial was developed using an example use case for a project being completed by the National Park Service. **The goal of the project is to use GEDI L1B Version 2 data to observe GEDI waveforms over Redwood National Park in northern California.** \n",
    "\n",
    "This tutorial will show how to use Python to open GEDI L1B Version 2 files, visualize the sub-orbit of GEDI points (shots), subset to a region of interest, visualize GEDI full waveforms, and export subsets of GEDI science dataset (SDS) layers as GeoJSON files that can be loaded into GIS and/or Remote Sensing software programs. \n",
    "- [Redwood National Park GeoJSON](https://github.com/nasa/GEDI-Data-Resources/Additional_files/RedwoodNP.geojson)  \n",
    "    - Contains the administrative boundary for Redwood National Park, available from: [Administrative Boundaries of National Park System Units 12/31/2017 - National Geospatial Data Asset (NGDA) NPS National Parks Dataset](https://irma.nps.gov/DataStore/DownloadFile/594958)\n",
    "\n",
    "***    \n",
    "### Data Used in the Example:  \n",
    "- **GEDI L1B Geolocated Waveform Data Global Footprint Level - [GEDI01_B.002](https://doi.org/10.5067/GEDI/GEDI01_B.002)**\n",
    "    - _The purpose of the L1B dataset is to provide geolocated waveforms and supporting datasets for each laser shot for all eight GEDI beams. This includes corrected and smoothed waveforms, geolocation parameters, and geophysical corrections._  \n",
    "     - **Science Dataset (SDS) layers:**  \n",
    "        - /geolocation/latitude_bin0  \n",
    "        - /geolocation/longitude_bin0  \n",
    "        - /shot_number  \n",
    "        - /stale_return_flag  \n",
    "        - /geolocation/degrade  \n",
    "        - /rx_sample_count  \n",
    "        - /rx_sample_start_index  \n",
    "        - /rxwaveform  \n",
    "        - /geolocation/elevation_bin0  \n",
    "        - /geolocation/elevation_lastbin  \n",
    "        - /geolocation/digital_elevation_model  \n",
    "\n",
    "***  \n",
    "# Topics Covered:\n",
    "1. [**Get Started**](#getstarted)  \n",
    "    1.1 Import Packages    \n",
    "    1.2 Set Up the Working Environment and Retrieve Files      \n",
    "2. [**Import and Interpret Data**](#importinterpret)      \n",
    "    2.1 Open a GEDI HDF5 File and Read File Metadata     \n",
    "    2.2 Read SDS Metadata and Subset by Beam   \n",
    "3. [**Visualize a GEDI Sub-Orbit**](#visualizeorbit)      \n",
    "    3.1 Subset by Layer and Create a Geodataframe   \n",
    "    3.2 Visualize a Geodataframe\n",
    "4. [**Subset and Visualize Waveforms**](#subsetviswaveforms)        \n",
    "    4.1 Import and Extract Waveforms   \n",
    "    4.2 Visualize Waveforms  \n",
    "5. [**Quality Filtering**](#quality)      \n",
    "6. [**Plot Profile Transects**](#plottransects)    \n",
    "    6.1 Subset Beam Transects    \n",
    "    6.2 Plot Waveform Transects   \n",
    "7. [**Spatial Visualization**](#spatialvisualization)    \n",
    "    7.1 Import, Subset, and Quality Filter All Beams      \n",
    "    7.2 Spatial Subsetting   \n",
    "    7.3 Visualize All Beams: Elevation  \n",
    "8. [**Export Subsets as GeoJSON Files**](#exportgeojson)               \n",
    "***\n",
    "# Before Starting this Tutorial:\n",
    "\n",
    "This tutorial requires a compatible Python Environment and GEDI L1B Version 2 granule from June 19, 2019 (orbit 02932, sub-orbit `02`) to download. To setup the Python environment and download the file, follow the steps in sections 1 and 2 of the [set-up instruction](https://github.com/nasa/GEDI-Data-Resources/Setup/setup_instructions.md).  \n",
    "***\n",
    "\n",
    "## Source Code used to Generate this Tutorial:\n",
    "The repository containing all of the required files is located at: https://github.com/nasa/GEDI-Data-Resources \n",
    "\n",
    "<div class=\"alert alert-block alert-warning\" >\n",
    "<b>NOTE:</b> This tutorial was developed for GEDI L1B Version 2 HDF-EOS5 files and should only be used for that product. </div>   "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# 1. Get Started <a id=\"getstarted\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.1 Import Packages <a id=\"1.1\"></a>\n",
    "#### Import the required packages and set the input/working directory to run this Jupyter Notebook locally."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import h5py\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import geopandas as gp\n",
    "import geoviews as gv\n",
    "from geoviews import opts, tile_sources as gvts\n",
    "import holoviews as hv\n",
    "gv.extension('bokeh', 'matplotlib')\n",
    "from shapely.geometry import Point\n",
    "import warnings\n",
    "from shapely.errors import ShapelyDeprecationWarning\n",
    "warnings.filterwarnings(\"ignore\", category=ShapelyDeprecationWarning) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.2 Set Up the Working Environment and Retrieve Files<a id=\"1.2\"></a>\n",
    "#### The input directory is defined as the current working directory. Note that you will need to have the Jupyter Notebook and example data (.h5 and .geojson) stored in this directory in order to execute the tutorial successfully."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "inDir = os.getcwd().split(\"python\")[0]   # Set input directory to the current working directory\n",
    "os.chdir(inDir)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### In this section, the GEDI .h5 file has been downloaded to the `inDir` defined above. You will need to download the file directly from the LP DAAC Data Pool in order to execute this tutorial.\n",
    "### Direct Link to file:\n",
    " - https://e4ftl01.cr.usgs.gov/GEDI/GEDI01_B.002/2019.06.19/GEDI01_B_2019170155833_O02932_02_T02267_02_005_01_V002.h5    \n",
    "   \n",
    "#### Make sure to download the file into the `inDir` directory defined above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gediFiles = [g for g in os.listdir(inDir+\"\\\\data\") if g.startswith('GEDI01_B') and g.endswith('.h5')]  # List all GEDI L1B .h5 files in inDir\n",
    "gediFiles"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# 2. Import and Interpret Data <a id=\"importinterpret\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.1 Open a GEDI HDF5 File and Read File Metadata <a id=\"2.1\"></a>\n",
    "#### Read the file using `h5py`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "L1B = f'data\\\\{gediFiles[0]}'\n",
    "L1B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The standard format  for GEDI Version 2 filenames is as follows:\n",
    "> **GEDI01_B**: Product Short Name    \n",
    "**2019170155833**: Julian Date and Time of Acquisition (YYYYDDDHHMMSS)  \n",
    "**O02932**: Orbit Number  \n",
    "**02**: Sub-Orbit Granule Number (1-4)  \n",
    "**T02267**: Track Number (Reference Ground Track)   \n",
    "**02**: Positioning and Pointing Determination System (PPDS) type (00 is predict, 01 rapid, 02 and higher is final)  \n",
    "**005**: PGE Version Number    \n",
    "**01**: Granule Production Version  \n",
    "**V002**: Product Version  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Read in a GEDI HDF5 file using the `h5py` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gediL1B = h5py.File(L1B, 'r')  # Read file using h5py"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Navigate the HDF5 file below. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(gediL1B.keys())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### The GEDI HDF5 file contains groups in which data and metadata are stored.\n",
    "#### First, the `METADATA` group contains the file-level metadata."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "list(gediL1B['METADATA'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This contains useful information such as the creation date, PGEVersion, and VersionID. Below, print the file-level metadata attributes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for g in gediL1B['METADATA']['DatasetIdentification'].attrs: print(g) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(gediL1B['METADATA']['DatasetIdentification'].attrs['purpose'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.2 Read SDS Metadata and Subset by Beam <a id=\"2.2\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "####  The GEDI instrument consists of 3 lasers producing a total of 8 beam ground transects. The eight remaining groups contain data for each of the eight GEDI beam transects. For additional information, be sure to check out: https://gedi.umd.edu/instrument/specifications/."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beamNames = [g for g in gediL1B.keys() if g.startswith('BEAM')]\n",
    "beamNames"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### One useful piece of metadata to retrieve from each beam transect is whether it is a full power beam or a coverage beam. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for g in gediL1B['BEAM0000'].attrs: print(g)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for b in beamNames: \n",
    "    print(f\"{b} is a {gediL1B[b].attrs['description']}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Below, pick one of the full power beams that will be used to retrieve GEDI L1B waveforms in Section 3. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beamNames = ['BEAM0110']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Identify all the objects in the GEDI HDF5 file below. \n",
    "Note: This step may take a while to complete."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gediL1B_objs = []\n",
    "gediL1B.visit(gediL1B_objs.append)                                           # Retrieve list of datasets\n",
    "gediSDS = [o for o in gediL1B_objs if isinstance(gediL1B[o], h5py.Dataset)]  # Search for relevant SDS inside data file\n",
    "[i for i in gediSDS if beamNames[0] in i][0:10]                              # Print the first 10 datasets for selected beam"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# 3. Visualize a GEDI Orbit <a id=\"visualizeorbit\"></a>\n",
    "#### In the section below, import GEDI L1B SDS layers into a `GeoPandas` GeoDataFrame for the beam specified above. \n",
    "#### Use the `latitude_bin0` and `longitude_bin0` to create a `shapely` point for each GEDI shot location. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1 Subset by Layer and Create a Geodataframe <a id=\"3.1\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Read in the SDS and take a representative sample (every 100th shot) and append to lists, then use the lists to generate a `pandas` dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "lonSample, latSample, shotSample, srfSample, degradeSample, beamSample = [], [], [], [], [], []  # Set up lists to store data\n",
    "\n",
    "# Open the SDS\n",
    "lats = gediL1B[f'{beamNames[0]}/geolocation/latitude_bin0'][()]\n",
    "lons = gediL1B[f'{beamNames[0]}/geolocation/longitude_bin0'][()]\n",
    "shots = gediL1B[f'{beamNames[0]}/shot_number'][()]\n",
    "srf = gediL1B[f'{beamNames[0]}/stale_return_flag'][()]\n",
    "degrade = gediL1B[f'{beamNames[0]}/geolocation/degrade'][()]\n",
    "\n",
    "# Take every 100th shot and append to list\n",
    "for i in range(len(shots)):\n",
    "    if i % 100 == 0:\n",
    "        shotSample.append(str(shots[i]))\n",
    "        lonSample.append(lons[i])\n",
    "        latSample.append(lats[i])\n",
    "        srfSample.append(srf[i])\n",
    "        degradeSample.append(degrade[i])\n",
    "        beamSample.append(beamNames[0])\n",
    "            \n",
    "# Write all of the sample shots to a dataframe\n",
    "latslons = pd.DataFrame({'Beam': beamSample, 'Shot Number': shotSample, 'Longitude': lonSample, 'Latitude': latSample,\n",
    "                         'Stale Return Flag': srfSample, 'Degrade': degradeSample})\n",
    "latslons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Above is a dataframe containing columns describing the beam, shot number, lat/lon location, and quality information about each shot."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Clean up variables that will no longer be needed\n",
    "del beamSample, degrade, degradeSample, gediL1B_objs, latSample, lats, lonSample, lons, shotSample, shots, srf, srfSample "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Below, create an additional column called 'geometry' that contains a `shapely` point generated from each lat/lon location from the shot. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Take the lat/lon dataframe and convert each lat/lon to a shapely point\n",
    "latslons['geometry'] = latslons.apply(lambda row: Point(row.Longitude, row.Latitude), axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Next, convert to a `Geopandas` GeoDataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert to a Geodataframe\n",
    "latslons = gp.GeoDataFrame(latslons)\n",
    "latslons = latslons.drop(columns=['Latitude','Longitude'])\n",
    "latslons['geometry']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Pull out and plot an example `shapely` point below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "latslons['geometry'][0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2 Visualize a GeoDataFrame <a id=\"3.2\"></a>\n",
    "#### In this section, use the GeoDataFrame and the `geoviews` python package to spatially visualize the location of the GEDI shots on a basemap and import a GeoJSON file of the spatial region of interest for the use case example: Redwood National Park."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define a function for visualizing GEDI points\n",
    "def pointVisual(features, vdims):\n",
    "    return (gvts.EsriImagery * gv.Points(features, vdims=vdims).options(tools=['hover'], height=500, width=900, size=5, \n",
    "                                                                        color='yellow', fontsize={'xticks': 10, 'yticks': 10, \n",
    "                                                                                                  'xlabel':16, 'ylabel': 16}))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Import a GeoJSON of Redwood National Park as an additional GeoDataFrame. Note that you will need to have downloaded the [GeoJSON](https://git.earthdata.nasa.gov/projects/LPDUR/repos/gedi-V2-tutorials/browse/RedwoodNP.geojson) from the bitbucket repo containing this tutorial and have it saved in the same directory as this Jupyter Notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "redwoodNP = gp.GeoDataFrame.from_file('Data\\\\RedwoodNP.geojson')  # Import GeoJSON as GeoDataFrame"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "redwoodNP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "redwoodNP['geometry'][0]  # Plot GeoDataFrame"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Defining the vdims below will allow you to hover over specific shots and view information about them."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a list of geodataframe columns to be included as attributes in the output map\n",
    "vdims = []\n",
    "for f in latslons:\n",
    "    if f not in ['geometry']:\n",
    "        vdims.append(f)\n",
    "vdims"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Below, combine a plot of the Redwood National Park Boundary (combine two `geoviews` plots using `*`) with the point visual mapping function defined above in order to plot (1) the representative GEDI shots, (2) the region of interest, and (3) a basemap layer. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Call the function for plotting the GEDI points\n",
    "gv.Polygons(redwoodNP['geometry']).opts(line_color='red', color=None) * pointVisual(latslons, vdims = vdims)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Above is a good illustration of the new GEDI _Version 2_ sub-orbit granules (remember that GEDI _Version 1_  files are stored as one ISS orbit). One of the benefits of using geoviews is the interactive nature of the output plots. Use the tools to the right of the map above to zoom in and find the shots intersecting Redwood National Park. \n",
    "> (**HINT**: find where the orbit intersects the west coast of the United States)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Below is a screenshot of the region of interest:\n",
    "![GEDI_L1B_V2_Tutorial_1.png](../../img/GEDI_L1B_V2_Tutorial_1.png \"Sample of GEDI L1B V2 shots in yellow (orbit 2932) plotted over Redwood National Park, USA.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Side Note: Wondering what the 0's and 1's for `stale_return_flag` and `degrade` mean?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"stale_return_flag: {gediL1B[b]['stale_return_flag'].attrs['description']}\")\n",
    "print(f\"degrade: {gediL1B[b]['geolocation']['degrade'].attrs['description']}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gediL1B[b]['stale_return_flag'].attrs.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"stale_return_flag: {gediL1B[b]['stale_return_flag'].attrs['description']}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### We will show an example of how to quality filter GEDI data in section 5.\n",
    "#### After finding one of the shots within Redwood NP, find the index for that shot number so that we can find the correct waveform to visualize in Section 4. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Each GEDI shot has a unique shot identifier (shot number) that is available within each data group of the product. The shot number is important to retain in any data subsetting as it will allow the user to link any shot record back to the original orbit data, and to link any shot and its data between the L1 and L2 products. The standard format  for GEDI Version 2 Shots is as follows:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Shot: 29320600200465601\n",
    "> **2932**: Orbit Number      \n",
    "**06**: Beam Number    \n",
    "**0**: Reserved for future use     \n",
    "**02**: Sub-orbit Granule Number    \n",
    "**004**: Minor frame number   \n",
    "**65601**: Shot index  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "shot = 29320600200465601"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "index = np.where(gediL1B[f'{beamNames[0]}/shot_number'][()]==shot)[0][0]  # Set the index for the shot identified above\n",
    "index"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del latslons  # No longer need the geodataframe used to visualize the full GEDI orbit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# 4. Subset and Visualize Waveforms <a id=\"subsetviswaveforms\"></a>\n",
    "#### In this section, learn how to extract and subset specific waveforms and plot them using `holoviews`. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.1 Import and Extract Waveforms<a id=\"4.1\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### In order to find and extract the full waveform for the exact index we are interested in, instead of importing the entire waveform dataset (over 100 million values!), we will use the `rx_sample_count` and `rx_sample_start_index` to identify the location of the waveform that we are interested in visualizing and just extract those values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# From the SDS list, use list comprehension to find sample_count, sample_start_index, and rxwaveform\n",
    "sdsCount = gediL1B[[g for g in gediSDS if g.endswith('/rx_sample_count') and beamNames[0] in g][0]]\n",
    "sdsStart = gediL1B[[g for g in gediSDS if g.endswith('/rx_sample_start_index') and beamNames[0] in g][0]]\n",
    "sdsWaveform = [g for g in gediSDS if g.endswith('/rxwaveform') and beamNames[0] in g][0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Print the description for each of these datasets to better understand how we will use them to extract specific waveforms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"rxwaveform is {gediL1B[sdsWaveform].attrs['description']}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Next, read how to use the rx_sample_count and rx_sample_start_index layers:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"rx_sample_count is {sdsCount.attrs['description']}\")\n",
    "print(f\"rx_sample_start_index is {sdsStart.attrs['description']}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Use `rx_sample_count` and `rx_sample_start_index` to identify the location of each waveform in `rxwaveform`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wfCount = sdsCount[index]           # Number of samples in the waveform\n",
    "wfStart = int(sdsStart[index] - 1)  # Subtract one because python array indexing starts at 0 not 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Next grab additional information about the shot, including the unique `shot_number`, and lat/lon location."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "wfShot = gediL1B[f'{beamNames[0]}/shot_number'][index]\n",
    "wfLat = gediL1B[f'{beamNames[0]}/geolocation/latitude_bin0'][index]\n",
    "wfLon = gediL1B[f'{beamNames[0]}/geolocation/longitude_bin0'][index]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Put everything together to identify the waveform we want to extract:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "print(f\"The waveform located at: {str(wfLat)}, {str(wfLon)} (shot ID: {wfShot}, index {index}) is from beam {beamNames[0]} \\\n",
    "      and is stored in rxwaveform beginning at index {wfStart} and ending at index {wfStart + wfCount}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### In order to plot a full waveform, you also need to import the elevation recorded at `bin0` (the first elevation recorded) and `lastbin` or the last elevation recorded for that waveform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Grab the elevation recorded at the start and end of the full waveform capture\n",
    "zStart = gediL1B[f'{beamNames[0]}/geolocation/elevation_bin0'][index]   # Height of the start of the rx window\n",
    "zEnd = gediL1B[f'{beamNames[0]}/geolocation/elevation_lastbin'][index]  # Height of the end of the rx window"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Extract the full waveform using the index start and count:\n",
    "#### Below you can see why it is important to extract the specific waveform that you are interested in: over 100 million values are stored in the rxwaveform dataset!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"{:,}\".format(gediL1B[sdsWaveform].shape[0]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Retrieve the waveform sds layer using the sample start index and sample count information to slice the correct dimensions\n",
    "waveform = gediL1B[sdsWaveform][wfStart: wfStart + wfCount]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.2 Visualize Waveforms\n",
    "#### Below, plot the extracted waveform using the elevation difference from bin0 to last_bin on the y axis, and the waveform energy returned or amplitude (digital number = dn) on the x axis. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find elevation difference from start to finish and divide into equal intervals based on sample_count\n",
    "zStretch = np.add(zEnd, np.multiply(range(wfCount, 0, -1), ((zStart - zEnd) / int(wfCount))))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# match the waveform amplitude values with the elevation and convert to Pandas df\n",
    "wvDF = pd.DataFrame({'Amplitude (DN)': waveform, 'Elevation (m)': zStretch})\n",
    "wvDF"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Next, use the `Holoviews` python package (hv) to plot the waveform."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hv.Curve(wvDF) # Basic line graph plotting the waveform"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Congratulations! You have plotted your first waveform.\n",
    "#### Above is a basic line plot showing amplitude (DN) as a function of elevation from the rxwaveform for the specific shot selected. Next, add additional chart elements and make the graph interactive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a holoviews interactive Curve plot with additional parameters defining the plot aesthetics \n",
    "wfVis = hv.Curve(wvDF).opts(color='darkgreen', tools=['hover'], height=500, width=400,\n",
    "           xlim=(np.min(waveform) - 10, np.max(waveform) + 10), ylim=(np.min(zStretch), np.max(zStretch)),\n",
    "           fontsize={'xticks':10, 'yticks':10,'xlabel':16, 'ylabel': 16, 'title':13}, line_width=2.5, title=f'{str(wfShot)}')\n",
    "wfVis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### As you can see, the selected shot does not look so interesting--this waveform shape looks more characteristic of some type of low canopy/bare ground than a tree canopy. If you zoom in to the GEDI shot it appears to be over the river itself or possibly a sand bar:\n",
    "![GEDI_L1B_V2_Tutorial_2.png](../../img/GEDI_L1B_V2_Tutorial_2.png \"GEDI L1B Shot over a river in Redwood National Park, USA.\")\n",
    "#### Next, plot a couple more waveforms and see if you can capture a shot over the forest."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latlons = {}                          # Set up a dictionary to hold multiple waveforms\n",
    "latlons[wfShot] = Point(wfLon,wfLat)  # Create shapely point and add to dictionary\n",
    "\n",
    "# Retain waveform and quality layers to be exported later\n",
    "latlonsWF = [waveform]\n",
    "latlonsEL = [zStretch]\n",
    "latlonsSRF = [gediL1B[f'{beamNames[0]}/stale_return_flag'][index]]\n",
    "latlonsD = [gediL1B[f'{beamNames[0]}/geolocation/degrade'][index]]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define which observation to examine by selecting the three indexes prior to the one used above\n",
    "index = 45729"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use rx_sample_count and rx_sample_start_index to identify the location of each waveform\n",
    "wfCount = sdsCount[index]                                             # Number of samples in the waveform\n",
    "wfStart = int(sdsStart[index] - 1)                                    # Subtract 1, index starts at 0 not 1\n",
    "wfShot = gediL1B[f'{beamNames[0]}/shot_number'][index]                # Unique Shot Number\n",
    "wfLat = gediL1B[f'{beamNames[0]}/geolocation/latitude_bin0'][index]   # Latitude\n",
    "wfLon = gediL1B[f'{beamNames[0]}/geolocation/longitude_bin0'][index]  # Longitude\n",
    "latlons[wfShot] = Point(wfLon,wfLat)                                  # Create shapely point and add to dictionary\n",
    "\n",
    "# Grab the elevation recorded at the start and end of the full waveform capture\n",
    "zStart = gediL1B[f'{beamNames[0]}/geolocation/elevation_bin0'][index]   # Height of the start of the rx window\n",
    "zEnd = gediL1B[f'{beamNames[0]}/geolocation/elevation_lastbin'][index]  # Height of the end of the rx window\n",
    "\n",
    "# Retrieve the waveform sds layer using the sample start index and sample count information to slice the correct dimensions\n",
    "waveform = gediL1B[sdsWaveform][wfStart: wfStart + wfCount]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find elevation difference from start to finish, divide into equal intervals\n",
    "zStretch = np.add(zEnd, np.multiply(range(wfCount, 0, -1), ((zStart - zEnd) / int(wfCount))))\n",
    "\n",
    "# match the waveform amplitude values with the elevation and convert to Pandas df\n",
    "wvDF = pd.DataFrame({'Amplitude (DN)': waveform, 'Elevation (m)': zStretch})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Append waveform and quality layers to be exported later\n",
    "latlonsWF.append(waveform)\n",
    "latlonsEL.append(zStretch)\n",
    "latlonsSRF.append(gediL1B[f'{beamNames[0]}/stale_return_flag'][index])\n",
    "latlonsD.append(gediL1B[f'{beamNames[0]}/geolocation/degrade'][index])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a holoviews interactive Curve plot with additional parameters defining the plot aesthetics \n",
    "visL1B1 = hv.Curve(wvDF).opts(color='darkgreen', tools=['hover'], height=500, width=400,\n",
    "           xlim=(np.min(waveform) - 10, np.max(waveform) + 10), ylim=(np.min(zStretch), np.max(zStretch)),\n",
    "           fontsize={'xticks':10, 'yticks':10,'xlabel':16, 'ylabel': 16, 'title':13}, line_width=2.5, title=f'{str(wfShot)}')\n",
    "visL1B1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Now that is starting to look more like a very tall, multi-layered tree canopy! Continue with the next shot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define which observation to examine\n",
    "index = 45730"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use rx_sample_count and rx_sample_start_index to identify the location of each waveform\n",
    "wfCount = sdsCount[index]                                             # Number of samples in the waveform\n",
    "wfStart = int(sdsStart[index] - 1)                                    # Subtract 1, index starts at 0 not 1\n",
    "wfShot = gediL1B[f'{beamNames[0]}/shot_number'][index]                # Unique Shot Number\n",
    "wfLat = gediL1B[f'{beamNames[0]}/geolocation/latitude_bin0'][index]   # Latitude\n",
    "wfLon = gediL1B[f'{beamNames[0]}/geolocation/longitude_bin0'][index]  # Longitude\n",
    "latlons[wfShot] = Point(wfLon,wfLat)                                  # Create shapely point and add to dictionary\n",
    "\n",
    "# Grab the elevation recorded at the start and end of the full waveform capture\n",
    "zStart = gediL1B[f'{beamNames[0]}/geolocation/elevation_bin0'][index]   # Height of the start of the rx window\n",
    "zEnd = gediL1B[f'{beamNames[0]}/geolocation/elevation_lastbin'][index]  # Height of the end of the rx window\n",
    "\n",
    "# Retrieve the waveform sds layer using the sample start index and sample count information to slice the correct dimensions\n",
    "waveform = gediL1B[sdsWaveform][wfStart: wfStart + wfCount]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find elevation difference from start to finish, divide into equal intervals\n",
    "zStretch = np.add(zEnd, np.multiply(range(wfCount, 0, -1), ((zStart - zEnd) / int(wfCount))))\n",
    "\n",
    "# match the waveform amplitude values with the elevation and convert to Pandas df\n",
    "wvDF = pd.DataFrame({'Amplitude (DN)': waveform, 'Elevation (m)': zStretch})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Append waveform and quality layers to be exported later\n",
    "latlonsWF.append(waveform)\n",
    "latlonsEL.append(zStretch)\n",
    "latlonsSRF.append(gediL1B[f'{beamNames[0]}/stale_return_flag'][index])\n",
    "latlonsD.append(gediL1B[f'{beamNames[0]}/geolocation/degrade'][index])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a holoviews interactive Curve plot with additional parameters defining the plot aesthetics \n",
    "visL1B2 = hv.Curve(wvDF).opts(color='darkgreen', tools=['hover'], height=500, width=400,\n",
    "           xlim=(np.min(waveform) - 10, np.max(waveform) + 10), ylim=(np.min(zStretch), np.max(zStretch)),\n",
    "           fontsize={'xticks':10, 'yticks':10,'xlabel':16, 'ylabel': 16, 'title':13}, line_width=2.5, title=f'{str(wfShot)}')\n",
    "visL1B2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define which observation to examine\n",
    "index = 45731"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use rx_sample_count and rx_sample_start_index to identify the location of each waveform\n",
    "wfCount = sdsCount[index]                                             # Number of samples in the waveform\n",
    "wfStart = int(sdsStart[index] - 1)                                    # Subtract 1, index starts at 0 not 1\n",
    "wfShot = gediL1B[f'{beamNames[0]}/shot_number'][index]                # Unique Shot Number\n",
    "wfLat = gediL1B[f'{beamNames[0]}/geolocation/latitude_bin0'][index]   # Latitude\n",
    "wfLon = gediL1B[f'{beamNames[0]}/geolocation/longitude_bin0'][index]  # Longitude\n",
    "latlons[wfShot] = Point(wfLon,wfLat)                                  # Create shapely point and add to dictionary\n",
    "\n",
    "# Grab the elevation recorded at the start and end of the full waveform capture\n",
    "zStart = gediL1B[f'{beamNames[0]}/geolocation/elevation_bin0'][index]   # Height of the start of the rx window\n",
    "zEnd = gediL1B[f'{beamNames[0]}/geolocation/elevation_lastbin'][index]  # Height of the end of the rx window\n",
    "\n",
    "# Retrieve the waveform sds layer using the sample start index and sample count information to slice the correct dimensions\n",
    "waveform = gediL1B[sdsWaveform][wfStart: wfStart + wfCount]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Find elevation difference from start to finish, divide into equal intervals\n",
    "zStretch = np.add(zEnd, np.multiply(range(wfCount, 0, -1), ((zStart - zEnd) / int(wfCount))))\n",
    "\n",
    "# match the waveform amplitude values with the elevation and convert to Pandas df\n",
    "wvDF = pd.DataFrame({'Amplitude (DN)': waveform, 'Elevation (m)': zStretch})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Append waveform and quality layers to be exported later\n",
    "latlonsWF.append(waveform)\n",
    "latlonsEL.append(zStretch)\n",
    "latlonsSRF.append(gediL1B[f'{beamNames[0]}/stale_return_flag'][index])\n",
    "latlonsD.append(gediL1B[f'{beamNames[0]}/geolocation/degrade'][index])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a holoviews interactive Curve plot with additional parameters defining the plot aesthetics \n",
    "visL1B3 = hv.Curve(wvDF).opts(color='darkgreen', tools=['hover'], height=500, width=400,\n",
    "           xlim=(np.min(waveform) - 10, np.max(waveform) + 10), ylim=(np.min(zStretch), np.max(zStretch)),\n",
    "           fontsize={'xticks':10, 'yticks':10,'xlabel':16, 'ylabel': 16, 'title':13}, line_width=2.5, title=f'{str(wfShot)}')\n",
    "visL1B3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Lastly, plot the transect of four consecutive waveforms:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# The \"+\" symbol will plot multiple saved holoviews plots together\n",
    "visL1B1.opts(width=240) + visL1B2.opts(width=240, labelled=[]) + visL1B3.opts(width=240, labelled=[]) + \\\n",
    "wfVis.opts(width=240, labelled=[])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Notice above moving west to east (left to right) along the transect you can see the elevation lowering and the tree canopy decreasing as it encounters the sandbar/river."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Below, use geoviews to plot the location of the four points to verify what is seen above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert dict to geodataframe\n",
    "latlons = gp.GeoDataFrame({'Shot Number': list(latlons.keys()),'rxwaveform Amplitude (DN)': latlonsWF,\n",
    "                           'rxwaveform Elevation (m)': latlonsEL, 'Stale Return Flag': latlonsSRF, 'Degrade': latlonsD,\n",
    "                           'geometry': list(latlons.values())})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latlons['Shot Number'] = latlons['Shot Number'].astype(str)  # Convert shot number from integer to string"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "vdims"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a list of geodataframe columns to be included as attributes in the output map\n",
    "vdims = []\n",
    "for f in latlons:\n",
    "    if f not in ['geometry']:\n",
    "        if 'rxwaveform' not in f:\n",
    "            vdims.append(f)\n",
    "\n",
    "# Plot the geodataframe\n",
    "gv.Polygons(redwoodNP['geometry']).opts(line_color='red', color=None) * pointVisual(latlons, vdims=vdims)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Above, zoom in on the yellow dots until you can get a clearer view of all four shots.\n",
    "#### In the screenshots below, we can match the location of each point with the associated waveform. It appears to confirm the hypothesis that the elevation is decreasing as the shots get closer to the river, and the forest canopy decreases until it is detecting the sand bar/river by the final waveform in the sequence.\n",
    "![GEDI_L1B_V2_Tutorial_3.png](../../img/GEDI_L1B_V2_Tutorial_3.png \"Plots of Four GEDI L1B Waveforms\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![GEDI_L1B_V2_Tutorial_4.png](../../img/GEDI_L1B_V2_Tutorial_4.png \"Map showing sequence of four GEDI L1B Waveforms in Redwood National Park, USA.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del wfVis, visL1B1, visL1B2, visL1B3, waveform, wvDF, zStretch, latlonsWF, latlonsEL, latlonsSRF, latlonsD"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# 5. Quality Filtering <a id=\"quality\"></a>\n",
    "#### Now that you have the desired layers imported as a dataframe for the selected shots, let's perform quality filtering.\n",
    "#### Below, remove any shots where the `stale_return_flag` and `degrade` flag is not set to 0.  \n",
    "#### Non-zero `stale_return_flag` indicates that a \"stale\" cue point from the coarse search algorithm is being used. The syntax of the line below finds the rows \"where\" the stale return flag is equal to 0 (better quality) in the dataframe. If a row (shot) does not meet the condition, set all values equal to `nan` for that row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "latlons = latlons.where(latlons['Stale Return Flag'] == 0)  # Set any stale returns to NaN"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Non-zero values for the `degrade` flag indicate the shot occured during a degraded period. Therefore in a same way, the rows \"where\" the `degrade` flag is equal to 0 (better quality) are found. If a row (shot) does not meet the condition, set all values equal to `nan` for that row."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latlons = latlons.where(latlons['Degrade'] == 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Below, drop all of the shots that did not pass the quality filtering standards outlined above from the `transectDF`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "latlons = latlons.dropna()  # Drop all of the rows (shots) that did not pass the quality filtering above"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(f\"Quality filtering complete, {len(latlons)} high quality shots remaining.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Good news! It looks like all four of the example waveforms passed the initial quality filtering tests. For additional information on quality filtering GEDI data, be sure to check out: https://lpdaac.usgs.gov/resources/faqs/#how-should-i-quality-filter-gedi-l1b-l2b-data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del latlons"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# 6. Plot Profile Transects <a id=\"plottransects\"></a>\n",
    "#### In this section, plot a transect subset using waveforms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.1 Subset Beam Transects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Subset down to a smaller transect centered on the waveforms analyzed in the sections above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(index)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Grab 50 points before and after the shot visualized above\n",
    "start = index - 50\n",
    "end = index + 50 \n",
    "transectIndex = np.arange(start, end, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6.2 Plot Waveform Transects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### In order to get an idea of the length of the beam transect that you are plotting, you can plot the x-axis as distance, which is calculated below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate along-track distance\n",
    "distance = np.arange(0.0, len(transectIndex) * 60, 60)  # GEDI Shots are spaced 60 m apart"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### In order to plot each vertical value for each waveform, you will need to reformat the data structure to match what is needed by `holoviews` Path() capabilities. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create a list of tuples containing Shot number, and each waveform value at height (z)\n",
    "wfList = []\n",
    "for s, i in enumerate(transectIndex):\n",
    "    # Basic quality filtering from section 5\n",
    "    if gediL1B['BEAM0110/geolocation/degrade'][i] == 0 and gediL1B['BEAM0110/stale_return_flag'][i] == 0:\n",
    "        zStart = gediL1B['BEAM0110/geolocation/elevation_bin0'][i]\n",
    "        zEnd = gediL1B['BEAM0110/geolocation/elevation_lastbin'][i]\n",
    "        zCount = sdsCount[i]\n",
    "        zStretch = np.add(zEnd, np.multiply(range(zCount, 0, -1), ((zStart - zEnd) / int(zCount))))\n",
    "        waveform = gediL1B[sdsWaveform][sdsStart[i]: sdsStart[i] + zCount]\n",
    "        waves = []\n",
    "        for z, w in enumerate(waveform):\n",
    "            waves.append((distance[s], zStretch[z],w))  # Append Distance (x), waveform elevation (y) and waveform amplitude (z)\n",
    "        wfList.append(waves)\n",
    "    else:\n",
    "        print(f\"Shot {s} did not pass quality filtering and will be excluded.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Good news again, it looks like all of the waveforms in our transect passed the quality filtering test.\n",
    "\n",
    "### Below, plot each waveform by using `holoviews` Path() function. This will plot each individual waveform value by distance, with the amplitude plotted in the third dimension in shades of green."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "# Plot the waveform values for the transect with the waveform amplitude being plotted as the 'c' or colormap\n",
    "l1bVis = hv.Path(wfList, vdims='Amplitude').options(color='Amplitude', clim=(200,405), cmap='Greens', line_width=20, width=950,\n",
    "                                                    height=500, clabel='Amplitude',colorbar=True, \n",
    "                                                    xlabel='Distance Along Transect (m)', ylabel='Elevation (m)', \n",
    "                                                    fontsize={'title':16, 'xlabel':16, 'ylabel': 16, 'xticks':12, 'yticks':12, \n",
    "                                                              'clabel':12, 'cticks':10}, colorbar_position='right',\n",
    "                                                    title='GEDI L1B Waveforms across Redwood National Park: June 19, 2019')\n",
    "l1bVis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### If you zoom in to a subset of the transect, you can get even better detail showing the structure of the canopy from the `rxwaveform` dataset, where denser portions of the canopy are seen in darker shades of green.\n",
    "![GEDI_L1B_V2_Tutorial_5.png](../../img/GEDI_L1B_V2_Tutorial_5.png \"Full waveform amplitude plotted for a transect covering Redwood National Park, USA.\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del waveform, waves, wfList, zStretch"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# 7. Spatial Visualization<a id=\"spatialvisualization\"></a>\n",
    "#### Section 7 combines many of the techniques learned above including how to import GEDI datasets, perform quality filtering, spatial subsetting, and visualization. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.1 Import, Subset, and Quality Filter All Beams"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Below, re-open the GEDI L1B observation--but this time, loop through and import data for all 8 of the GEDI beams."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beamNames = [g for g in gediL1B.keys() if g.startswith('BEAM')]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beamNames"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Loop through each of the desired datasets (SDS) for each beam, append to lists, and transform into a `pandas` DataFrame."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set up lists to store data\n",
    "shotNum, dem, zElevation, zHigh, zLat, zLon, srf, degrade, beamI, rxCount, rxStart = ([] for i in range(11))  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Loop through each beam and open the SDS needed\n",
    "for b in beamNames:\n",
    "    [shotNum.append(h) for h in gediL1B[[g for g in gediSDS if g.endswith('/shot_number') and b in g][0]][()]]\n",
    "    [dem.append(h) for h in gediL1B[[g for g in gediSDS if g.endswith('/digital_elevation_model') and b in g][0]][()]]\n",
    "    [zElevation.append(h) for h in gediL1B[[g for g in gediSDS if g.endswith('/elevation_bin0') and b in g][0]][()]]  \n",
    "    [zHigh.append(h) for h in gediL1B[[g for g in gediSDS if g.endswith('/elevation_lastbin') and b in g][0]][()]]  \n",
    "    [zLat.append(h) for h in gediL1B[[g for g in gediSDS if g.endswith('/latitude_bin0') and b in g][0]][()]]  \n",
    "    [zLon.append(h) for h in gediL1B[[g for g in gediSDS if g.endswith('/longitude_bin0') and b in g][0]][()]]  \n",
    "    [srf.append(h) for h in gediL1B[[g for g in gediSDS if g.endswith('/stale_return_flag') and b in g][0]][()]]  \n",
    "    [degrade.append(h) for h in gediL1B[[g for g in gediSDS if g.endswith('/degrade') and b in g][0]][()]]  \n",
    "    [rxCount.append(h) for h in gediL1B[[g for g in gediSDS if g.endswith('/rx_sample_count') and b in g][0]][()]]  \n",
    "    [rxStart.append(h) for h in gediL1B[[g for g in gediSDS if g.endswith('/rx_sample_start_index') and b in g][0]][()]]  \n",
    "    [beamI.append(h) for h in [b] * len(gediL1B[[g for g in gediSDS if g.endswith('/shot_number') and b in g][0]][()])]  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Notice that we did not import the rxwaveform dataset--due to the large size of that dataset, we will wait until after spatial subsetting to import the waveforms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert lists to Pandas dataframe\n",
    "allDF = pd.DataFrame({'Shot Number': shotNum, 'Beam': beamI, 'Latitude': zLat, 'Longitude': zLon, 'Tandem-X DEM': dem,\n",
    "                      'Elevation bin 0 (m)': zElevation, 'Elevation last bin (m)': zHigh, 'Stale Return Flag': srf,\n",
    "                      'Degrade Flag': degrade, 'rx Start Index': rxStart, 'rx Sample Count': rxCount})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del shotNum, dem, zElevation, zHigh, zLat, zLon, srf, degrade, beamI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.2 Spatial Subsetting\n",
    "#### Below, subset the pandas dataframe using a simple bounding box region of interest. If you are interested in spatially clipping GEDI shots to a GeoJSON region of interest, be sure to check out the GEDI-Subsetter python script available at: https://git.earthdata.nasa.gov/projects/LPDUR/repos/gedi-subsetter/browse."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(allDF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Almost 800,000 shots are contained in this single GEDI sub-orbit granule! Below subset down to only the shots falling within this small bounding box encompassing Redwood National Park. `RedwoodNP` our `geopandas` geodataframe can be called for the \"envelope\" or smallest bounding box encompassing the entire region of interest. Here, use that as the bounding box for subsetting the GEDI shots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "redwoodNP.envelope[0].bounds"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "minLon, minLat, maxLon, maxLat = redwoodNP.envelope[0].bounds  # Define the min/max lat/lon from the bounds of Redwood NP"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Filter by the bounding box, which is done similarly to filtering by quality in section 5 above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "allDF = allDF.where(allDF['Latitude'] > minLat)\n",
    "allDF = allDF.where(allDF['Latitude'] < maxLat)\n",
    "allDF = allDF.where(allDF['Longitude'] > minLon)\n",
    "allDF = allDF.where(allDF['Longitude'] < maxLon)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### NOTE: It may take up to a few minutes to run the cell above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "allDF = allDF.dropna()  # Drop shots outside of the ROI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "len(allDF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Notice you have drastically reduced the number of shots you are working with (which will greatly enhance your experience in plotting them below). But first, remove any poor quality shots that exist within the ROI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Set any poor quality returns to NaN\n",
    "allDF = allDF.where(allDF['Stale Return Flag'].ne(1))\n",
    "allDF = allDF.where(allDF['Degrade Flag'] < 1)\n",
    "allDF = allDF.dropna()\n",
    "len(allDF)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Next create a `Shapely` Point out of each shot and insert it as the geometry column in the [soon to be geo]dataframe."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Take the lat/lon dataframe and convert each lat/lon to a shapely point\n",
    "allDF['geometry'] = allDF.apply(lambda row: Point(row.Longitude, row.Latitude), axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Convert to geodataframe\n",
    "allDF = gp.GeoDataFrame(allDF)\n",
    "allDF = allDF.drop(columns=['Latitude','Longitude'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7.3 Visualize All Beams: Elevation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Now, using the `pointVisual` function defined in section 3.2, plot the `geopandas` GeoDataFrame using `geoviews`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "allDF['Shot Number'] = allDF['Shot Number'].astype(str)  # Convert shot number to string\n",
    "\n",
    "vdims = []\n",
    "for f in allDF:\n",
    "    if f not in ['geometry']:\n",
    "        vdims.append(f)\n",
    "\n",
    "visual = pointVisual(allDF, vdims = vdims)\n",
    "visual * gv.Polygons(redwoodNP['geometry']).opts(line_color='red', color=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Feel free to pan and zoom in to the GEDI shots in yellow. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Now let's not only plot the points in the geodataframe but also add a colormap for Elevation bin0 (m). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "(gvts.EsriImagery * gv.Points(allDF, vdims=vdims).options(color='Elevation bin 0 (m)',cmap='terrain', size=3, tools=['hover'],\n",
    "                                                          clim=(min(allDF['Elevation bin 0 (m)']), \n",
    "                                                          max(allDF['Elevation bin 0 (m)'])), colorbar=True, clabel='Meters',\n",
    "                                                          title='GEDI Elevation bin0 over Redwood National Park: June 19, 2019',\n",
    "                                                          fontsize={'xticks': 10, 'yticks': 10, 'xlabel':16, 'clabel':12,\n",
    "                                                                    'cticks':10,'title':16,'ylabel':16})).options(height=500,\n",
    "                                                                                                                  width=900)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---\n",
    "# 8. Export Subsets as GeoJSON Files<a id=\"exportgeojson\"></a>\n",
    "#### In this section, export the GeoDataFrame as a `.geojson` file that can be easily opened in your favorite remote sensing and/or GIS software and will include an attribute table with all of the shots/values for each of the SDS layers in the dataframe."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### First, import the waveforms for the subset of shots to be exported."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "waves = []  # List to store each list of waveform values\n",
    "\n",
    "# Loop through each row (shot) in the data frame, extract waveform values, append to list of lists\n",
    "for _, shot in allDF.iterrows():\n",
    "    starts = int(shot['rx Start Index'])\n",
    "    wave = gediL1B[f\"{shot['Beam']}/rxwaveform\"][starts: starts + int(shot['rx Sample Count'])]\n",
    "    \n",
    "    # In order to export as GeoJSON, the type needs to be converted from a numpy array to a single comma separated string\n",
    "    waves.append(','.join([str(q) for q in wave]))\n",
    "    \n",
    "# Add to data frame\n",
    "allDF['rx Waveform'] = waves"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gediL1B.filename  # L1B Filename"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "outName = gediL1B.filename.replace('.h5', '.json')  # Create an output file name using the input file name\n",
    "outName"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "allDF.to_file(outName, driver='GeoJSON')  # Export to GeoJSON"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "del allDF, waves"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Success! You have now learned how to start working with GEDI L1B Version 2 files in Python as well as some interesting strategies for visualizing those data in order to better understand your specific region of interest. Using this Jupyter Notebook as a workflow, you should now be able to switch to GEDI files over your specific region of interest and re-run the notebook. Good Luck!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<div class=\"alert alert-block alert-info\">\n",
    "    <h1> Contact Information </h1>\n",
    "    <h3> Material written by LP DAAC<sup>1</sup> </h3>\n",
    "    <ul>\n",
    "        <b>Contact:</b> LPDAAC@usgs.gov <br> \n",
    "        <b>Voice:</b> +1-605-594-6116 <br>\n",
    "        <b>Organization:</b> Land Processes Distributed Active Archive Center (LP DAAC) <br>\n",
    "        <b>Website:</b> https://lpdaac.usgs.gov/ <br>\n",
    "        <b>Date last modified:</b> 03-27-2023 <br>\n",
    "    </ul>\n",
    "    \n",
    "<sup>1</sup>KBR Inc., contractor to the U.S. Geological Survey, Earth Resources Observation and Science (EROS) Center, Sioux Falls, South Dakota, 57198-001, USA. Work performed under USGS contract G15PD00467 for LP DAAC<sup>2</sup>.\n",
    "\n",
    "<sup>2</sup>LP DAAC Work performed under NASA contract NNG14HH33I.\n",
    "</div>"
   ]
  }
 ],
 "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.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}