Getting Started with GEDI L2A Data in Python

Details:

Published: May 13, 2020

⚠ Attention: The Global Ecosystem Dynamics Investigation (GEDI) Version 1 data products are no longer available for distribution from NASA’s Earthdata Search or LP DAAC’s Data Pool. Users are advised to use the most recent version of GEDI.

Getting Started with GEDI L2A Data in Python

This tutorial demonstrates how to work with the Elevation and Height Metrics (GEDI02_A.001) data product.

The Global Ecosystem Dynamics Investigation (GEDI) 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 products. The L1B and L2 GEDI products are archived and distributed in the HDF-EOS5 file format.


Use Case Example:

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 L2A data to observe tree canopy height and profile over Redwood National Park in northern California.

This tutorial will show how to use Python to open GEDI L2A files, visualize the full orbit of GEDI points (shots), subset to a region of interest, visualize GEDI canopy height, and export subsets of GEDI science dataset (SDS) layers as GeoJSON files that can be loaded into GIS and/or Remote Sensing software programs.


Data Used in the Example:

  • GEDI L2A Elevation and Height Metrics Data Global Footprint Level - GEDI02_A.001
    • The purpose of the L2A dataset is to provide waveform interpretation and extracted products from each GEDI waveform. This includes ground elevation, canopy top height, relative return energy metrics (describing canopy vertical structure, for example), and many other interpreted products from the return waveforms.
    • Science Dataset (SDS) layers:
      • /rh
      • /elev_lowestmode
      • /elev_highestreturn
      • /lat_lowestmode
      • /lon_lowestmode
      • /shot_number
      • /quality_flag
      • /digital_elevation_model
      • /degrade_flag
      • /sensitivity

Topics Covered:

  1. Get Started
    1.1 Import Packages
    1.2 Set Up the Working Environment and Retrieve Files
  2. Import and Interpret Data
    2.1 Open a GEDI HDF5 File and Read File Metadata
    2.2 Read SDS Metadata and Subset by Beam
  3. Visualize a GEDI Orbit
    3.1 Subset by Layer and Create a Geodataframe
    3.2 Visualize a Geodataframe
  4. Work with GEDI L2A Data
    4.1 Import and Extract Specific Shots
    4.2 Plot Relative Height Metrics
    4.3 Combine RH Metrics and Waveforms
    4.4 Select Data from non-Default Algorithm
  5. Plot Transects
    5.1 Quality Filtering
    5.2 Plot Beam Transects 5.3 Subset Beam Transects
    5.4 Plot RH Metrics Transects
  6. Spatial Visualization
    6.1 Import, Subset, and Quality Filter all Beams
    6.2 Spatial Subsetting
    6.3 Visualize All Beams: Canopy Height and Elevation
  7. Export Subsets as GeoJSON Files

Before Starting this Tutorial:

Setup and Dependencies

It is recommended to use Conda, an environment manager to set up a compatible Python environment. Download Conda for your OS here: https://www.anaconda.com/download/. Once you have Conda installed, Follow the instructions below to successfully setup a Python environment on Linux, MacOS, or Windows.

This Python Jupyter Notebook tutorial has been tested using Python version 3.7. Conda was used to create the python environment.

  • Using your preferred command line interface (command prompt, terminal, cmder, etc.) type the following to successfully create a compatible python environment:

    conda create -n geditutorial -c conda-forge --yes python=3.7 h5py shapely geopandas pandas geoviews holoviews

conda activate geditutorial

jupyter notebook

If you do not have Jupyter Notebook installed, you may need to run:

conda install jupyter notebook

Having trouble getting a compatible Python environment set up? Contact LP DAAC User Services.

If you prefer to not install Conda, the same setup and dependencies can be achieved by using another package manager such as pip.


Example Data:

This tutorial uses the GEDI L2A observation from June 19, 2019 (orbit 02932). Use the link below to download the file directly from the LP DAAC Data Pool:
- https://e4ftl01.cr.usgs.gov/GEDI/GEDI02_A.001/2019.06.19/GEDI02_A_2019170155833_O02932_T02267_02_001_01.h5 (5.93 GB)

A NASA Earthdata Login account is required to download the data used in this tutorial. You can create an account at the link provided.

You will need to have the file above downloaded into the same directory as this Jupyter Notebook in order to successfully run the code below.

Source Code used to Generate this Tutorial:

The repository containing all of the required files is located at: https://git.earthdata.nasa.gov/projects/LPDUR/repos/gedi-tutorials/browse
- Jupyter Notebook
- Redwood National Park GeoJSON
- 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
- waveform.csv

NOTE: This tutorial was developed for GEDI L2A files and should only be used for those products.

1. Get Started

1.1 Import Packages

Import the required packages and set the input/working directory to run this Jupyter Notebook locally.

import os
import h5py
import numpy as np
import pandas as pd
import geopandas as gp
from shapely.geometry import Point
import geoviews as gv
from geoviews import opts, tile_sources as gvts
import holoviews as hv
gv.extension('bokeh', 'matplotlib')

1.2 Set Up the Working Environment and Retrieve Files

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.

inDir = os.getcwd() + os.sep  # Set input directory to the current working directory
NOTE: If you have downloaded the tutorial materials to a different directory than the Jupyter Notebook, `inDir` above needs to be changed. You will also need to add a line: `os.chdir(inDir)` and execute it below.

In this section, a GEDI L2A .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.

Direct Link to file:

  • https://e4ftl01.cr.usgs.gov/GEDI/GEDI02_A.001/2019.06.19/GEDI02_A_2019170155833_O02932_T02267_02_001_01.h5

Make sure to download the file into the inDir directory defined above.

gediFiles = [g for g in os.listdir() if g.startswith('GEDI02_A') and g.endswith('.h5')]  # List GEDI L2A .h5 files in the inDir
gediFiles
['GEDI02_A_2019170155833_O02932_T02267_02_001_01.h5']

2. Import and Interpret Data

In this section, begin working with the GEDI Level 2A Elevation and Height Metrics Data. Begin by checking out the example L2A file metadata.

2.1 Open a GEDI HDF5 File and Read File Metadata

Read the file using h5py.

L2A = 'GEDI02_A_2019170155833_O02932_T02267_02_001_01.h5'
L2A
'GEDI02_A_2019170155833_O02932_T02267_02_001_01.h5'

The standard format for GEDI filenames is as follows:

GEDI02_A: Product Short Name
2019170155833: Julian Date and Time of Acquisition (YYYYDDDHHMMSS)
O02932: Orbit Number
T02267: Track Number
02: Positioning and Pointing Determination System (PPDS) type (00 is predict, 01 rapid, 02 and higher is final)
001: GOC SDS (software) release number
01: Granule Production Version

Read in a GEDI HDF5 file using the h5py package.

gediL2A = h5py.File(L2A, 'r')  # Read file using h5py

Navigate the HDF5 file below.

list(gediL2A.keys())
['BEAM0000',
 'BEAM0001',
 'BEAM0010',
 'BEAM0011',
 'BEAM0101',
 'BEAM0110',
 'BEAM1000',
 'BEAM1011',
 'METADATA']

The GEDI HDF5 file contains groups in which data and metadata are stored.

First, the METADATA group contains the file-level metadata.

list(gediL2A['METADATA'])
['DatasetIdentification']

This contains useful information such as the creation date, PGEVersion, and VersionID. Below, print the file-level metadata attributes.

for g in gediL2A['METADATA']['DatasetIdentification'].attrs: print(g)
PGEVersion
VersionID
abstract
characterSet
creationDate
credit
fileName
language
originatorOrganizationName
purpose
shortName
spatialRepresentationType
status
topicCategory
uuid
print(gediL2A['METADATA']['DatasetIdentification'].attrs['purpose'])
The purpose of the L2A dataset is to provide waveform interpretation and extracted products from each GEDI waveform.  This includes ground elevation, canopy top height, relative return energy metrics (describing canopy vertical structure, for example), and many other interpreted products from the return waveforms.

2.2 Read SDS Metadata and Subset by Beam

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/.

beamNames = [g for g in gediL2A.keys() if g.startswith('BEAM')]
beamNames
['BEAM0000',
 'BEAM0001',
 'BEAM0010',
 'BEAM0011',
 'BEAM0101',
 'BEAM0110',
 'BEAM1000',
 'BEAM1011']

One useful piece of metadata to retrieve from each beam transect is whether it is a full power beam or a coverage beam.

for g in gediL2A['BEAM0000'].attrs: print(g)
description
for b in beamNames:
    print(f"{b} is a {gediL2A[b].attrs['description']}")
BEAM0000 is a Coverage beam
BEAM0001 is a Coverage beam
BEAM0010 is a Coverage beam
BEAM0011 is a Coverage beam
BEAM0101 is a Full power beam
BEAM0110 is a Full power beam
BEAM1000 is a Full power beam
BEAM1011 is a Full power beam

Below, pick one of the full power beams that will be used to retrieve GEDI L2A relative height metrics in Section 3.

beamNames = ['BEAM0110']

Identify all the objects in the GEDI HDF5 file below.

Note: This step may take a while to complete.

gediL2A_objs = []
gediL2A.visit(gediL2A_objs.append)                                           # Retrieve list of datasets
gediSDS = [o for o in gediL2A_objs if isinstance(gediL2A[o], h5py.Dataset)]  # Search for relevant SDS inside data file
[i for i in gediSDS if beamNames[0] in i][0:10]                              # Print the first 10 datasets for selected beam
['BEAM0110/ancillary/l2a_alg_count',
 'BEAM0110/beam',
 'BEAM0110/channel',
 'BEAM0110/degrade_flag',
 'BEAM0110/delta_time',
 'BEAM0110/digital_elevation_model',
 'BEAM0110/elev_highestreturn',
 'BEAM0110/elev_lowestmode',
 'BEAM0110/elevation_bias_flag',
 'BEAM0110/elevation_bin0_error']

3. Visualize a GEDI Orbit

In the section below, import GEDI L2A SDS layers into a GeoPandas GeoDataFrame for the beam specified above.

Use the lat_lowestmode and lon_lowestmode to create a shapely point for each GEDI shot location.

3.1 Subset by Layer and Create a Geodataframe

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.

lonSample, latSample, shotSample, qualitySample, beamSample = [], [], [], [], []  # Set up lists to store data

# Open the SDS
lats = gediL2A[f'{beamNames[0]}/lat_lowestmode'][()]
lons = gediL2A[f'{beamNames[0]}/lon_lowestmode'][()]
shots = gediL2A[f'{beamNames[0]}/shot_number'][()]
quality = gediL2A[f'{beamNames[0]}/quality_flag'][()]

# Take every 100th shot and append to list
for i in range(len(shots)):
    if i % 100 == 0:
        shotSample.append(str(shots[i]))
        lonSample.append(lons[i])
        latSample.append(lats[i])
        qualitySample.append(quality[i])
        beamSample.append(beamNames[0])

# Write all of the sample shots to a dataframe
latslons = pd.DataFrame({'Beam': beamSample, 'Shot Number': shotSample, 'Longitude': lonSample, 'Latitude': latSample,
                         'Quality Flag': qualitySample})
latslons
Beam Shot Number Longitude Latitude Quality Flag
0 BEAM0110 29320618800000001 111.996300 -51.803868 0
1 BEAM0110 29320604600000101 112.039132 -51.803905 0
2 BEAM0110 29320614600000201 112.080271 -51.803836 0
3 BEAM0110 29320600400000301 112.121445 -51.803737 0
4 BEAM0110 29320610400000401 112.162622 -51.803621 0
... ... ... ... ... ...
9792 BEAM0110 29320617400979201 88.208452 -51.803578 0
9793 BEAM0110 29320603200979301 88.249610 -51.803614 0
9794 BEAM0110 29320613200979401 88.290753 -51.803581 0
9795 BEAM0110 29320623200979501 88.331913 -51.803548 0
9796 BEAM0110 29320609000979601 88.373089 -51.803506 0

9797 rows × 5 columns

Above is a dataframe containing columns describing the beam, shot number, lat/lon location, and quality information about each shot.

# Clean up variables that will no longer be needed
del beamSample, quality, qualitySample, gediL2A_objs, latSample, lats, lonSample, lons, shotSample, shots

Below, create an additional column called 'geometry' that contains a shapely point generated from each lat/lon location from the shot.

# Take the lat/lon dataframe and convert each lat/lon to a shapely point
latslons['geometry'] = latslons.apply(lambda row: Point(row.Longitude, row.Latitude), axis=1)

Next, convert to a Geopandas GeoDataFrame.

# Convert to a Geodataframe
latslons = gp.GeoDataFrame(latslons)
latslons = latslons.drop(columns=['Latitude','Longitude'])
latslons['geometry']
0       POINT (111.99630 -51.80387)
1       POINT (112.03913 -51.80391)
2       POINT (112.08027 -51.80384)
3       POINT (112.12145 -51.80374)
4       POINT (112.16262 -51.80362)
                   ...             
9792     POINT (88.20845 -51.80358)
9793     POINT (88.24961 -51.80361)
9794     POINT (88.29075 -51.80358)
9795     POINT (88.33191 -51.80355)
9796     POINT (88.37309 -51.80351)
Name: geometry, Length: 9797, dtype: geometry

Pull out and plot an example shapely point below.

latslons['geometry'][0]

Single Shapely point.

3.2 Visualize a GeoDataFrame

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.

# Define a function for visualizing GEDI points
def pointVisual(features, vdims):
    return (gvts.EsriImagery * gv.Points(features, vdims=vdims).options(tools=['hover'], height=500, width=900, size=5,
                                                                        color='yellow', fontsize={'xticks': 10, 'yticks': 10,
                                                                                                  'xlabel':16, 'ylabel': 16}))

Import a geojson of Redwood National Park as an additional GeoDataFrame. Note that you will need to have downloaded the geojson from the bitbucket repo containing this tutorial and have it saved in the same directory as this Jupyter Notebook.

redwoodNP = gp.GeoDataFrame.from_file('RedwoodNP.geojson')  # Import geojson as GeoDataFrame
redwoodNP
GIS_LOC_ID UNIT_CODE GROUP_CODE UNIT_NAME UNIT_TYPE META_MIDF LANDS_CODE DATE_EDIT GIS_NOTES geometry
0 None REDW None Redwood National Park None None None Shifted 0.06 miles MULTIPOLYGON (((-124.01829 41.44539, -124.0184...
redwoodNP['geometry'][0]  # Plot GeoDataFrame

Shapefile of Redwood National Park

Defining the vdims below will allow you to hover over specific shots and view information about them.

# Create a list of geodataframe columns to be included as attributes in the output map
vdims = []
for f in latslons:
    if f not in ['geometry']:
        vdims.append(f)
vdims
['Beam', 'Shot Number', 'Quality Flag']

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.

# Call the function for plotting the GEDI points
gv.Polygons(redwoodNP['geometry']).opts(line_color='red', color=None) * pointVisual(latslons, vdims = vdims)

Global GEDI Orbit Map.

Above is a good illustration of the full GEDI orbit (GEDI 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.

(HINT: find where the orbit intersects the west coast of the United States)

Below is a screenshot of the region of interest:

alt text

Side Note: Wondering what the 0's and 1's for quality_flag mean?

print(f"Quality Flag: {gediL2A[b]['quality_flag'].attrs['description']}")
Quality Flag: Flag simpilfying selection of most useful data

Above, 0 is poor quality and a quality_flag value of 1 indicates the laser shot meets criteria based on energy, sensitivity, amplitude, and real-time surface tracking quality. We will show an example of how to quality filter GEDI data in section 6.1.

After finding one of the shots within Redwood NP, find the index for that shot number so that we can find the correct RH metrics to visualize in Section 4.

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 Shots is as follows:

Shot: 29320619900465601

2932: Orbit Number
06: Beam Number
199: Minor frame number (0-241)
00465601: Shot number within orbit

del latslons  # No longer need the geodataframe used to visualize the full GEDI orbit

4. Work with GEDI L2A Data

Interpretation of RH Metrics

The GEDI L2A data product provides relative height (RH) metrics, which are “lidar perceived” metrics that have the following characteristics:

1. RH100 = elev_highestreturn - elev_lowestmode
2. The RH metrics are intended for vegetated surfaces. Results over bare/water surfaces are still valid but may present some confusing results.
3. The lower RH metrics (e.g., RH10) will often have negative values, particularly in low canopy cover conditions. This is because a relatively high fraction of the waveform energy is from the ground and below elev_lowestmode. For example, if the ground return contains 30% of the energy, then RH1 through 15 are likely to be below 0 since half of the ground energy from the ground return is below the center of the ground return, which is used to determine the mean ground elevation in the footprint (elev_lowestmode).

In this section, import and extract a specific GEDI L2A shot. From there, plot relative height (RH) metrics, and then combine with the L1B full waveform data. Note that we will continue to use BEAM0110 from Section 3.

4.1 Import and Extract Specific Shots

Notice that there are thousands of datasets available in the GEDI L2A product. In the code blocks below, you will subset to just a few of the datasets available.

len(gediSDS)
4232
beamNames
['BEAM0110']

Begin by subsetting just to the selected high power beam:

beamSDS = [g for g in gediSDS if beamNames[0] in g]  # Subset to a single beam
len(beamSDS)
529

We will set the shot index used as an example from the GEDI L1B Tutorial to show how to subset a single shot of GEDI L2A data.

shot = 29320619500465599
index = np.where(gediL2A[f'{beamNames[0]}/shot_number'][()]==shot)[0][0]  # Set the index for the shot identified above
index
465598

4.2 Plot Relative Height Metrics

In section 4.2, import the relative height metrics and begin exploring how to plot them.

rh = gediL2A[[g for g in beamSDS if g.endswith('/rh')][0]]  # Relative Height Metrics

Print the description for the rh dataset to better understand relative height (RH) metrics.

print(f"rh is {rh.attrs['description']}")
rh is Relative height metrics at 1 % interval

Next, bring in other useful L2A datasets such as lat_lowestmode and lon_lowestmode.

# Bring in the desired SDS
lats = gediL2A[f'{beamNames[0]}/lat_lowestmode'][()]  # Latitude
lons = gediL2A[f'{beamNames[0]}/lon_lowestmode'][()]  # Longitude

Grab the location and relative height metrics for the shot defined above:

rhLat = lats[index]
rhLon = lons[index]
rhShot1 = rh[index]

Put everything together to identify the shot that we want to extract:

print(f"The shot is located at: {str(rhLat)}, {str(rhLon)} (shot ID: {shot}, index {index}) and is from beam {beamNames[0]}.")
The shot is located at: 41.28472739326018, -124.03109998658007 (shot ID: 29320619500465599, index 465598) and is from beam BEAM0110.

In order to plot the RH metrics, you also need to import the elevation recorded at elev_lowestmode (the first elevation recorded) and elev_highestreturn or the last elevation recorded for that waveform.

# Grab the elevation recorded at the start and end of the RH metrics
zElevation = gediL2A[[g for g in beamSDS if g.endswith('/elev_lowestmode')][0]][index]  # Elevation
zTop = gediL2A[[g for g in beamSDS if g.endswith('/elev_highestreturn')][0]][index]     # Elevation Highest Return

Below, convert canopy height to canopy elevation by adding the ground elevation to each RH 1% interval.

rhShot = [z + zElevation for z in rhShot1]  # To convert canopy height to canopy elevation, add the elevation to each value
rh25 = rhShot[24]                           # 25%
rh50 = rhShot[49]                           # 50%  
rh75 = rhShot[74]                           # 75%

Next, use the holoviews (hv) package to make a line plot showing the elevation at each percentile for the selected shot.

rhVis = hv.Curve(rhShot, label='Default Algorithm (a1)')
rhVis = rhVis.opts(color='black', tools=['hover'], height=500, width=400, title='GEDI L2A Relative Height Metrics',
                   xlabel='Percent Energy Returned', ylabel='Elevation (m)', xlim=(0,100),ylim=(np.min(rhShot),np.max(rhShot)),
                   fontsize={'title':14, 'xlabel':16, 'ylabel': 16, 'legend': 14, 'xticks':12, 'yticks':12}, line_width=3.5)
rhVis

Relative cumulative RH profile.

Congratulations! You have plotted your first relative cumulative RH profile.

The rh dataset stores the relative height metrics at 1% intervals, and thus each shot contains 101 values representing rh at 0-100%.

Next, add some additional context to the cumulative RH profile, including metrics such as the elev_lowestmode (ground elevation), the elev_highestreturn (highest reflecting surface height), and the relative height metrics at each quartile.

# Create plots for L2A Metrics
zX = [0,100]                   # set up list from 0 to 100 to create the line
zY = [zElevation, zElevation]  # ground elevation
zT = [zTop, zTop]              # highest return

# Set up plots for each of the desired values
zVis = hv.Curve((zX, zY), label='Ground Return').opts(color='saddlebrown', tools=['hover'], height=550, width=400, line_width=2)
ztVis = hv.Curve((zX, zT), label='RH100').opts(color='navy', tools=['hover'], height=550, width=400, line_width=2)
rh25Vis = hv.Curve((zX, [rh25,rh25]),label='RH25').opts(color='lightblue',tools=['hover'], height=550, width=400, line_width=2)
rh50Vis = hv.Curve((zX, [rh50,rh50]),label='RH50').opts(color='mediumblue',tools=['hover'], height=550, width=400, line_width=2)
rh75Vis = hv.Curve((zX, [rh75,rh75]),label='RH75').opts(color='darkblue',tools=['hover'], height=550, width=400, line_width=2)
# Plot all of the metrics together
l2aVis = rhVis * zVis * ztVis * rh25Vis * rh50Vis * rh75Vis
l2aVis.opts(show_legend=True, legend_position='bottom_right', title='GEDI L2A Relative Height Metrics', ylabel='Elevation (m)',
                   xlabel='Percent Energy Returned', xlim=(0, 100), ylim=(np.min(rhShot) + 1.5, np.max(rhShot) + 5), height=600,
                   fontsize={'title':16, 'xlabel':16, 'ylabel': 16, 'legend': 14, 'xticks':12, 'yticks':12}, width=400)

Graphic of GEDI L2A Relative Height Metrics

Based on the figure above, it looks like the densest portion of the tree canopy for this shot is at around 45-55 meters.

4.3 Combine RH Metrics and Waveforms

Note that you will need to have downloaded the waveform.csv file from the repo in order to execute the section below.

In this section, we will import the corresponding waveform for our selected shot.

wvDF = pd.read_csv('waveform.csv')
wvDF
Amplitude (DN) Elevation (m)
0 227.00992 111.252897
1 226.57410 111.103165
2 226.66390 110.953434
3 227.26736 110.803702
4 228.11644 110.653971
... ... ...
1241 230.74078 -74.564040
1242 230.04378 -74.713771
1243 229.24507 -74.863503
1244 228.71117 -75.013235
1245 228.53528 -75.162966

1246 rows × 2 columns

Next, create a holoviews curve plot of the waveform. If you are interested in learning more about plotting L1B waveforms, be sure to check out the GEDI L1B Tutorial.

# Create a holoviews interactive Curve plot with additional parameters defining the plot aesthetics
visL1B = hv.Curve(wvDF).opts(color='darkgreen', tools=['hover'], height=600, width=400,
           xlim=(np.min(wvDF['Amplitude (DN)']) - 10, np.max(wvDF['Amplitude (DN)']) + 10),
           ylim=(np.min(wvDF['Elevation (m)']), np.max(wvDF['Elevation (m)'])),
           fontsize={'xticks':10, 'yticks':10,'xlabel':16, 'ylabel': 16, 'title':13}, line_width=2.5, title=f'{str(shot)}')
visL1B

Curve plot of the waveform.

Next, combine the L2A Relative Height (RH) metrics to the corresponding L1B waveform.

visL1B.opts(height=600, width=400, ylim=(np.min(rhShot), np.max(rhShot)+5), ylabel='Elevation (m)', xlabel='Amplitude (DN)') \
+ l2aVis.opts(height=600, width=400, ylim=(np.min(rhShot), np.max(rhShot)+5))

combined  L2A Relative Height (RH) metrics to the corresponding L1B waveform.

Above, notice the amplitude spike where the ground return is defined, and the close grouping of RH50 and RH75 over the dense portion of the upper canopy detected for this shot.

4.4 Select Data from the Non-Default Algorithm

The GEDI L2A product provides waveform processing results for multiple algorithm settings. Results for a default algorithm selection are provided in the root directory of the data product for each beam. This is currently set to the output of algorithm setting group 1 (see Table 5 of the GEDI L02 User Guide) and will be updated as post-launch cal/val progresses. Algorithm Setting Group 1 is the algorithm used for the data that we extracted in the sections above.

Elevation and height metrics outputs for all algorithm setting groups can be found in the geolocation subgroup of the L2A data product. For example:
- elev_lowestreturn_a is the elevation of lowest return detected using algorithm setting group - rh_a are the relative height metrics at 1% intervals using algorithm (in cm).

In some cases, the selection of an alternative algorithm setting will provide a better result. See the GEDI L02A User Guide for additional information.

geolocationSDS = [b for b in beamSDS if '/geolocation/' in b]  # Select all datasets within the geolocation subgroup

Next, select the datasets for Algorithm Setting Group 2 (a2).

a2SDS = [g for g in geolocationSDS if g.endswith('a2')]  # Select algorithm 2 datasets
a2SDS
['BEAM0110/geolocation/elev_highestreturn_a2',
 'BEAM0110/geolocation/elev_lowestmode_a2',
 'BEAM0110/geolocation/elev_lowestreturn_a2',
 'BEAM0110/geolocation/elevs_allmodes_a2',
 'BEAM0110/geolocation/energy_lowestmode_a2',
 'BEAM0110/geolocation/lat_highestreturn_a2',
 'BEAM0110/geolocation/lat_lowestmode_a2',
 'BEAM0110/geolocation/lat_lowestreturn_a2',
 'BEAM0110/geolocation/lats_allmodes_a2',
 'BEAM0110/geolocation/lon_highestreturn_a2',
 'BEAM0110/geolocation/lon_lowestmode_a2',
 'BEAM0110/geolocation/lon_lowestreturn_a2',
 'BEAM0110/geolocation/lons_allmodes_a2',
 'BEAM0110/geolocation/num_detectedmodes_a2',
 'BEAM0110/geolocation/quality_flag_a2',
 'BEAM0110/geolocation/rh_a2',
 'BEAM0110/geolocation/sensitivity_a2']

Repeat Section 4.2 for Algorithm Setting Group 2.

# Bring in the desired SDS
lata2 = gediL2A[[a for a in a2SDS if a.endswith('/lat_lowestmode_a2')][0]][index]          # Latitude
lona2 = gediL2A[[a for a in a2SDS if a.endswith('/lon_lowestmode_a2')][0]][index]          # Latitude
rhShot1 = gediL2A[[a for a in a2SDS if a.endswith('/rh_a2')][0]][index] / 100              # Relative height metrics
zElevationa2 = gediL2A[[a for a in a2SDS if a.endswith('/elev_lowestmode_a2')][0]][index]  # Elevation
zTopa2 = gediL2A[[a for a in a2SDS if a.endswith('/elev_highestreturn_a2')][0]][index]     # Elevation Highest Return
rhShota2 = [z + zElevationa2 for z in rhShot1]  # To convert canopy height to canopy elevation, add the elevation to each value

Notice above that the non-default algorithm RH metrics are stored in cm, so we convert to meters by dividing by 100.

rhVisA2 = hv.Curve(rhShota2, label='Algorithm Setting Group 2 (a2)')
rhVisA2 = rhVisA2.opts(color='green', tools=['hover'], height=500, width=400, title='GEDI L2A Relative Height Metrics (A2)',
                   xlabel='Percent Energy Returned', ylabel='Elevation (m)', xlim=(0,100),ylim=(np.min(rhShot),np.max(rhShot)),
                   fontsize={'title':14, 'xlabel':16, 'ylabel': 16, 'legend': 14, 'xticks':12, 'yticks':12}, line_width=3.5)
rhVisA2

GEDI L2A Relative Height Metrics (A2)

Next, combine with the relative height metrics from Algorithm Setting Group 1:

(rhVis * rhVisA2).opts(show_legend=True, legend_position='bottom_right', title='GEDI L2A Relative Height Metrics by Algorithm',
                       ylabel='Elevation (m)', xlabel='Percent Energy Returned', xlim=(0, 100),
                       ylim=(np.min(rhShot) - 1.5, np.max(rhShot) + 1.5), height=600, width=600,
                   fontsize={'title':16, 'xlabel':16, 'ylabel': 16, 'legend': 14, 'xticks':12, 'yticks':12})

GEDI L2A Relative Height Metrics by Algorithm

Above we can see slight differences in the RH values between the two algorithms, particularly in the lower percentiles.


5. Plot Transects

Next, import a number of desired SDS layers for BEAM0110 (for the entire orbit) and create a pandas Dataframe to store the arrays.

# Open all of the desired SDS
dem = gediL2A[[g for g in beamSDS if g.endswith('/digital_elevation_model')][0]][()]
zElevation = gediL2A[[g for g in beamSDS if g.endswith('/elev_lowestmode')][0]][()]
zHigh = gediL2A[[g for g in beamSDS if g.endswith('/elev_highestreturn')][0]][()]
zLat = gediL2A[[g for g in beamSDS if g.endswith('/lat_lowestmode')][0]][()]
zLon = gediL2A[[g for g in beamSDS if g.endswith('/lon_lowestmode')][0]][()]
rh = gediL2A[[g for g in beamSDS if g.endswith('/rh')][0]][()]
quality = gediL2A[[g for g in beamSDS if g.endswith('/quality_flag')][0]][()]
degrade = gediL2A[[g for g in beamSDS if g.endswith('/degrade_flag')][0]][()]
sensitivity = gediL2A[[g for g in beamSDS if g.endswith('/sensitivity')][0]][()]
shotNums = gediL2A[f'{beamNames[0]}/shot_number'][()]

# Create a shot index
shotIndex = np.arange(shotNums.size)
canopyHeight = [r[100] for r in rh]  # Grab RH100 (index 100 for each RH metrics)
# Take the DEM, GEDI-produced Elevation, and RH Metrics and add to a Pandas dataframe
transectDF = pd.DataFrame({'Shot Index': shotIndex, 'Shot Number': shotNums, 'Latitude': zLat, 'Longitude': zLon,
                           'Tandem-X DEM': dem, 'Elevation (m)': zElevation, 'Canopy Elevation (m)': zHigh,
                           'Canopy Height (rh100)': canopyHeight, 'Quality Flag': quality, 'Degrade Flag': degrade,
                           'Sensitivity': sensitivity})
transectDF
Shot Index Shot Number Latitude Longitude Tandem-X DEM Elevation (m) Canopy Elevation (m) Canopy Height (rh100) Quality Flag Degrade Flag Sensitivity
0 0 29320618800000001 -51.803868 111.996300 -999999.0 21242.515625 21242.515625 0.0 0 0 -3.436965
1 1 29320618900000002 -51.803867 111.996712 -999999.0 21242.505859 21242.505859 0.0 0 0 30.496670
2 2 29320619000000003 -51.803867 111.997123 -999999.0 21242.496094 21242.496094 0.0 0 0 8.071431
3 3 29320619100000004 -51.803867 111.997535 -999999.0 21242.484375 21242.484375 0.0 0 0 -212.896439
4 4 29320619200000005 -51.803866 111.997946 -999999.0 21242.474609 21242.474609 0.0 0 0 -6.853874
... ... ... ... ... ... ... ... ... ... ... ...
979694 979694 29320618400979695 -51.803445 88.411747 -999999.0 18017.906250 18017.906250 0.0 0 0 22.138037
979695 979695 29320618500979696 -51.803445 88.412159 -999999.0 18017.296875 18017.296875 0.0 0 0 4.475757
979696 979696 29320618600979697 -51.803444 88.412570 -999999.0 18017.884766 18017.884766 0.0 0 0 10.112548
979697 979697 29320618700979698 -51.803444 88.412981 -999999.0 18017.275391 18017.275391 0.0 0 0 424.691803
979698 979698 29320618800979699 -51.803443 88.413393 -999999.0 18017.263672 18017.263672 0.0 0 0 15.887813

979699 rows × 11 columns

Notice the unusual values listed above--those shots are flagged as poor quality and will be removed in Section 5.1.

Now that you have the desired SDS into a pandas dataframe, begin plotting the entire beam transect:

# Plot Canopy Height
canopyVis = hv.Scatter((transectDF['Shot Index'], transectDF['Canopy Height (rh100)']))
canopyVis.opts(color='darkgreen', height=500, width=900, title=f'GEDI L2A Full Transect {beamNames[0]}',
               fontsize={'title':16, 'xlabel':16, 'ylabel': 16}, size=0.1, xlabel='Shot Index', ylabel='Canopy Height (m)')

GEDI L2A Full Transect BEAM0110

Congratulations! You have plotted your first GEDI full orbit beam transect. Notice above that things look a little messy--before we dive deeper into plotting full transects, let's quality filter the shots in the section below.

del canopyVis, canopyHeight, degrade, dem, quality, sensitivity, shotIndex, shotNums, zElevation, zHigh, zLat, zLon

5.1 Quality Filtering

Now that you have the desired layers imported as a dataframe for the entire beam transect, let's perform quality filtering.

Below, remove any shots where the quality_flag is set to 0 by defining those shots as nan.

The syntax of the line below can be read as: in the dataframe, find the rows "where" the quality flag is not equal (ne) to 0. If a row (shot) does not meet the condition, set all values equal to nan for that row.

transectDF = transectDF.where(transectDF['Quality Flag'].ne(0))  # Set any poor quality returns to NaN
transectDF
Shot Index Shot Number Latitude Longitude Tandem-X DEM Elevation (m) Canopy Elevation (m) Canopy Height (rh100) Quality Flag Degrade Flag Sensitivity
0 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
1 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
2 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
3 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
4 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
... ... ... ... ... ... ... ... ... ... ... ...
979694 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
979695 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
979696 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
979697 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
979698 NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

979699 rows × 11 columns

Below, quality filter even further by using the degrade_flag (Greater than zero if the shot occurs during a degrade period, zero otherwise) and the Sensitivity layer, using a threshold of 0.95.

From the GEDI Level 2 User Guide:

Use the sensitivity metric available in L2A and L2B to select "best" data. The L2A and L2B quality_flag datasets use a conservative sensitivity threshold of 0.9 over land (0.5 over ocean), but under some conditions (e.g. dense forest) the user may benefit from selecting a higher threshold.

transectDF = transectDF.where(transectDF['Degrade Flag'].ne(1))
transectDF = transectDF.where(transectDF['Sensitivity'] > 0.95)

Below, drop all of the shots that did not pass the quality filtering standards outlined above from the transectDF.

transectDF = transectDF.dropna()  # Drop all of the rows (shots) that did not pass the quality filtering above
print(f"Quality filtering complete, {len(transectDF)} high quality shots remaining.")
Quality filtering complete, 79429 high quality shots remaining.

5.2 Plot Beam Transects

Next, plot the full remaining transect of high quality values using holoviews Scatter(). Combine the Tandem-X derived elevation, the GEDI-derived elevation, and the Canopy Top Elevation in a combined holoviews plot.

# Remove any shots where there is no Tandem-X Elevation values
transectDF = transectDF.where(transectDF['Tandem-X DEM'] != -999999.0).dropna()

# Plot Digital Elevation Model
demVis = hv.Scatter((transectDF['Shot Index'], transectDF['Tandem-X DEM']), label='Tandem-X DEM')
demVis = demVis.opts(color='black', height=500, width=900, fontsize={'xlabel':16, 'ylabel': 16}, size=1.5)
# Plot GEDI-Retrieved Elevation
zVis = hv.Scatter((transectDF['Shot Index'], transectDF['Elevation (m)']), label='GEDI-derived Elevation')
zVis = zVis.opts(color='saddlebrown', height=500, width=900, fontsize={'xlabel':16, 'ylabel': 16}, size=1.5)
# Plot Canopy Top Elevation
rhVis = hv.Scatter((transectDF['Shot Index'], transectDF['Canopy Elevation (m)']), label='Canopy Top Elevation')
rhVis = rhVis.opts(color='darkgreen', height=500, width=900, fontsize={'xlabel':16, 'ylabel': 16}, size=1.5,
                   tools=['hover'], xlabel='Shot Index', ylabel='Elevation (m)')
# Combine all three scatterplots
(demVis * zVis * rhVis).opts(show_legend=True, legend_position='top_left',fontsize={'title':15, 'xlabel':16, 'ylabel': 16},
                             title=f'{beamNames[0]} Full Transect: {L2A.split(".")[0]}')

BEAM0110 Full Transect

The plot still looks a bit messy this far zoomed out--feel free to pan, zoom, and explore different areas of the plot. The shot plotted in section 4 was at index 46598. If you zoom into the high-quality shots between 4.000e+5 and 5.000e+5, you will find the portion of the transect intersecting Redwood National Park, seen below:

alt text

5.3 Subset Beam Transects

Now, subset down to a smaller transect centered on the GEDI shot analyzed in the sections above.

print(index)
465598
# Grab 50 points before and after the shot visualized above
start = index - 50
end = index + 50
print(f"The transect begins at ({transectDF['Latitude'][start]}, {transectDF['Longitude'][start]}) and ends at ({transectDF['Latitude'][end]}, {transectDF['Longitude'][end]}).")
The transect begins at (41.26951477815523, -124.05868759659765) and ends at (41.299873598763384, -124.00358737366548).

Below, subset the transect using .loc.

transectDF = transectDF.loc[start:end]  # Subset the Dataframe to only the selected region of interest over Redwood NP

5.4 Plot RH Metrics Transects

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.

# Calculate along-track distance
distance = np.arange(0.0, len(transectDF.index) * 60, 60)  # GEDI Shots are spaced 60 m apart
transectDF['Distance'] = distance                          # Add Distance as a new column in the dataframe

But how do you determine where the ground is? Next, create a similar plot but using the rh metrics from the L2A file. Notice that the rh metrics contain 100 values representing each percentile, so below we iterate through and append each value with the correct distance and elevation value to the list.

# Create lists of canopy height, distance, and canopy elevation
rhList, distList, ceList = [], [], []
for s, i in enumerate(transectDF['Shot Index']):
    i = int(i)
    rhi = rh[i]
    for r in rhi:
        rhList.append(r)
        distList.append(distance[s])
        ceList.append(r + transectDF['Elevation (m)'][i])  # Convert to Canopy Top Elevation

Now plot canopy height vs. shot index as a scatter plot.

title = 'Relative Height Metrics (0-100) across Redwood National Park: June 19, 2019'
hv.Scatter((distList,rhList)).opts(color='darkgreen', alpha=0.1, height=500, width=900, size=7.5,
                                                xlabel='Distance Along Transect (m)', ylabel='Canopy Height (m)', title=title,
                                                fontsize={'title':16, 'xlabel':16, 'ylabel': 16, 'xticks':12, 'yticks':12})

Relative Height Metrics (0-100) across Redwood National Park: June 19, 2019

Above, it looks like the transect is covering some tall trees! However, a quick google search confirms that Redwood trees can top out at over 100 m! The graph above is showing canopy height, but now add elevation to get a better idea of the overall three-dimensional structure of this transect across Redwood National Park.

# Plot canopy elevation instead of canopy height
rhVis = hv.Scatter((distList,ceList), label='RH 0-100').opts(color='darkgreen', alpha=0.1, height=500, width=900, size=5)

Next, add in the ground elevation.

# Plot GEDI-Retrieved Elevation
demVis = hv.Curve((transectDF['Distance'], transectDF['Elevation (m)']), label='Ground Elevation')
demVis = demVis.opts(color='black', height=500, width=900, line_width=3, xlabel='Distance Along Transect (m)',
                     ylabel='Canopy Height (m)', fontsize={'title':14, 'xlabel':16, 'ylabel': 16, 'xticks':12, 'yticks':12},
                     title='Relative Height Metrics and Elevation across Redwood National Park: June 19, 2019')

# Plot Elevation and RH metrics together
(demVis * rhVis).opts(legend_position='bottom_right')

Relative Height Metrics and Elevation across Redwood National Park: June 19, 2019

Above, you can get an idea about the terrain over the region of interest, particularly the classic "V" representing the river valley that is bisected by the transect. In terms of vegetation structure, this plot does a good job of showing not only which portions of the canopy are taller, but also where they are denser (darker shades of green).

del ceList, distList, rhList, rhi
del beamSDS, distance, rh, transectDF

At this point you have visualized the elevation, canopy, and vertical structure of specific footprints over Redwood National Park, and for a transect cutting through the national park. In section 6 you will look at mapping all of the high-quality shots from all eight GEDI beams for a given region of interest in order to gain knowledge on the spatial distribution and characteristics of the canopy over Redwood National Park.


6. Spatial Visualization

Section 6 combines many of the techniques learned above including how to import GEDI datasets, perform quality filtering, spatial subsetting, and visualization.

6.1 Import, Subset, and Quality Filter All Beams

Below, re-open the GEDI L2A observation--but this time, loop through and import data for all 8 of the GEDI beams.

beamNames = [g for g in gediL2A.keys() if g.startswith('BEAM')]
beamNames
['BEAM0000',
 'BEAM0001',
 'BEAM0010',
 'BEAM0011',
 'BEAM0101',
 'BEAM0110',
 'BEAM1000',
 'BEAM1011']

Loop through each of the desired datasets (SDS) for each beam, append to lists, and transform into a pandas DataFrame.

# Set up lists to store data
shotNum, dem, zElevation, zHigh, zLat, zLon, rh25, rh98, rh100 ,quality ,degrade, sensitivity ,beamI = ([] for i in range(13))  

Below, note that we are going to pull out the relative height metrics at the 25th, 98th, and 100th (canopy height) percentiles.

# Loop through each beam and open the SDS needed
for b in beamNames:
    [shotNum.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/shot_number') and b in g][0]][()]]
    [dem.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/digital_elevation_model') and b in g][0]][()]]
    [zElevation.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/elev_lowestmode') and b in g][0]][()]]  
    [zHigh.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/elev_highestreturn') and b in g][0]][()]]  
    [zLat.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/lat_lowestmode') and b in g][0]][()]]  
    [zLon.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/lon_lowestmode') and b in g][0]][()]]  
    [rh25.append(h[25]) for h in gediL2A[[g for g in gediSDS if g.endswith('/rh') and b in g][0]][()]]  
    [rh98.append(h[98]) for h in gediL2A[[g for g in gediSDS if g.endswith('/rh') and b in g][0]][()]]
    [rh100.append(h[100]) for h in gediL2A[[g for g in gediSDS if g.endswith('/rh') and b in g][0]][()]]  
    [quality.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/quality_flag') and b in g][0]][()]]  
    [degrade.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/degrade_flag') and b in g][0]][()]]  
    [sensitivity.append(h) for h in gediL2A[[g for g in gediSDS if g.endswith('/sensitivity') and b in g][0]][()]]  
    [beamI.append(h) for h in [b] * len(gediL2A[[g for g in gediSDS if g.endswith('/shot_number') and b in g][0]][()])]  

Note that the cell above may take a few seconds to complete.

# Convert lists to Pandas dataframe
allDF = pd.DataFrame({'Shot Number': shotNum, 'Beam': beamI, 'Latitude': zLat, 'Longitude': zLon, 'Tandem-X DEM': dem,
                      'Elevation (m)': zElevation, 'Canopy Elevation (m)': zHigh, 'Canopy Height (rh100)': rh100, 'RH 98': rh98,
                      'RH 25': rh25, 'Quality Flag': quality, 'Degrade Flag': degrade, 'Sensitivity': sensitivity})
del beamI, degrade, dem, gediSDS, rh100, rh98, rh25, quality, sensitivity, zElevation, zHigh, zLat, zLon, shotNum

6.2 Spatial Subsetting

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.

len(allDF)
3547051

Over 3.5 million shots are contained in this single GEDI orbit! 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.

redwoodNP.envelope[0].bounds
(-124.16015705494489,
 41.080601363502545,
 -123.84950230520286,
 41.83981133687605)
minLon, minLat, maxLon, maxLat = redwoodNP.envelope[0].bounds  # Define the min/max lat/lon from the bounds of Redwood NP

Filter by the bounding box, which is done similarly to filtering by quality in section 6.1 above.

allDF = allDF.where(allDF['Latitude'] > minLat)
allDF = allDF.where(allDF['Latitude'] < maxLat)
allDF = allDF.where(allDF['Longitude'] > minLon)
allDF = allDF.where(allDF['Longitude'] < maxLon)

NOTE: It may take up to a few minutes to run the cell above.

allDF = allDF.dropna()  # Drop shots outside of the ROI
len(allDF)
4477

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.

# Set any poor quality returns to NaN
allDF = allDF.where(allDF['Quality Flag'].ne(0))
allDF = allDF.where(allDF['Degrade Flag'].ne(1))
allDF = allDF.where(allDF['Sensitivity'] > 0.95)
allDF = allDF.dropna()
len(allDF)
2081

Down to roughly 2000 shots, next create a Shapely Point out of each shot and insert it as the geometry column in the [soon to be geo]dataframe.

# Take the lat/lon dataframe and convert each lat/lon to a shapely point
allDF['geometry'] = allDF.apply(lambda row: Point(row.Longitude, row.Latitude), axis=1)
# Convert to geodataframe
allDF = gp.GeoDataFrame(allDF)
allDF = allDF.drop(columns=['Latitude','Longitude'])

6.3 Visualize All Beams: Canopy Height and Elevation

Now, using the pointVisual function defined in section 3.2, plot the geopandas GeoDataFrame using geoviews.

allDF['Shot Number'] = allDF['Shot Number'].astype(str)  # Convert shot number to string

vdims = []
for f in allDF:
    if f not in ['geometry']:
        vdims.append(f)

visual = pointVisual(allDF, vdims = vdims)
visual * gv.Polygons(redwoodNP['geometry']).opts(line_color='red', color=None)

GEDI shots in yellow.

Feel free to pan and zoom in to the GEDI shots in yellow.

Now let's plot the points in the geodataframe and add a colormap for Canopy Height (m) and Elevation (m).

# Plot the basemap and geoviews Points, defining the color as the Canopy Height for each shot
(gvts.EsriImagery * gv.Points(allDF, vdims=vdims).options(color='Canopy Height (rh100)',cmap='plasma', size=3, tools=['hover'],
                                                          clim=(0,102), colorbar=True, clabel='Meters',
                                                          title='GEDI Canopy Height over Redwood National Park: June 19, 2019',
                                                          fontsize={'xticks': 10, 'yticks': 10, 'xlabel':16, 'clabel':12,
                                                                    'cticks':10,'title':16,'ylabel':16})).options(height=500,
                                                                                                                  width=900)

GEDI Canopy Height over Redwood National Park: June 19, 2019

Above and in the screenshot below, notice the higher canopy heights (shades of yellow) over the Redwood stands of the national park vs. other types of forests (pink-blue) vs. the low-lying (and consequently flat) profiles over lakes and rivers (purple).

alt text

Next, take a look at the GEDI-derived elevation over the shots. Notice below that the colormap is changed to 'terrain'.

(gvts.EsriImagery * gv.Points(allDF, vdims=vdims).options(color='Elevation (m)',cmap='terrain', size=3, tools=['hover'],
                                                          clim=(min(allDF['Elevation (m)']), max(allDF['Elevation (m)'])),
                                                          colorbar=True, clabel='Meters',
                                                          title='GEDI Elevation over Redwood National Park: June 19, 2019',
                                                          fontsize={'xticks': 10, 'yticks': 10, 'xlabel':16, 'clabel':12,
                                                                    'cticks':10,'title':16,'ylabel':16})).options(height=500,
                                                                                                                  width=900)

GEDI Elevation over Redwood National Park: June 19, 2019

Success! You have now learned how to start working with GEDI L2A 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!


7. Export Subsets as GeoJSON Files

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.

gediL2A.filename  # L2A Filename
'GEDI02_A_2019170155833_O02932_T02267_02_001_01.h5'
outName = gediL2A.filename.replace('.h5', '.json')  # Create an output file name using the input file name
outName
'GEDI02_A_2019170155833_O02932_T02267_02_001_01.json'
allDF.to_file(outName, driver='GeoJSON')  # Export to GeoJSON
del allDF

Contact Information

Material written by Cole Krehbiel1

    Contact: LPDAAC@usgs.gov
    Voice: +1-605-594-6116
    Organization: Land Processes Distributed Active Archive Center (LP DAAC)
    Website: https://lpdaac.usgs.gov/
    Date last modified: 05-11-2021
1KBR 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 DAAC2. 2LP DAAC Work performed under NASA contract NNG14HH33I.

Relevant Products

Product Long Name
GEDI01_B.001 GEDI L1B Geolocated Waveform Data Global Footprint Level
GEDI02_A.001 GEDI L2A Elevation and Height Metrics Data Global Footprint Level

Tools

Name Filters Description