Getting Started with GEDI L2B Data in Python¶
This tutorial demonstrates how to work with the Canopy Cover and Vertical Profile Metrics (GEDI02_B.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$^{o}$ N and 51.6$^{o}$ 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 L2B data to observe tree canopy height, cover, and profile over Redwood National Park in northern California.
This tutorial will show how to use Python to open GEDI L2B files, visualize the full orbit of GEDI points (shots), subset to a region of interest, visualize GEDI canopy height and vertical profile metrics, 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 L2B Canopy Cover and Vertical Profile Metrics Data Global Footprint Level - GEDI02_B.001
- The purpose of the L2B dataset is to extract biophysical metrics from each GEDI waveform. These metrics are based on the directional gap probability profile derived from the L1B waveform and include canopy cover, Plant Area Index (PAI), Plant Area Volume Density (PAVD) and Foliage Height Diversity (FHD).
- Science Dataset (SDS) layers:
- /geolocation/digital_elevation_model
- /geolocation/elev_lowestmode
- /geolocation/elev_highestreturn
- /geolocation/lat_lowestmode
- /geolocation/lon_lowestmode
- /rh100
- /l2b_quality_flag
- /degrade_flag
- /sensitivity
- /pai
- /pavd_z
- /geolocation/shot_number
- /dz
Topics Covered:¶
- Get Started
1.1 Import Packages
1.2 Set Up the Working Environment and Retrieve Files - Import and Interpret Data
2.1 Open a GEDI HDF5 File and Read File Metadata
2.2 Read SDS Metadata and Subset by Beam - Visualize a GEDI Orbit
3.1 Subset by Layer and Create a Geodataframe
3.2 Visualize a Geodataframe - Work with GEDI L2B Data
4.1 Import and Extract PAVD
4.2 Visualize PAVD - Work with GEDI L2B Beam Transects
5.1 Quality Filtering
5.2 Plot Beam Transects
5.3 Subset Beam Transects - Plot Profile Transects
6.1 Plot PAVD Transects - Spatial Visualization
7.1 Import, Subset, and Quality Filter all Beams
7.2 Spatial Subsetting
7.3 Visualize All Beams: Canopy Height, Elevation, and PAI - 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 holoviewsconda activate geditutorialjupyter 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 at: https://lpdaac.usgs.gov/lpdaac-contact-us/¶
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 L2B observation from June 19, 2019 (orbit 02932). Use the links below to download the files directly from the LP DAAC Data Pool:
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
NOTE: This tutorial was developed for GEDI L2B HDF-EOS5 files and should only be used for that product.
1. Get Started ¶
In [1]:
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')
In [2]:
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 .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_B.001/2019.06.19/GEDI02_B_2019170155833_O02932_T02267_02_001_01.h5
#### Make sure to download the files into theinDirdirectory defined above.
In [3]:
gediFiles = [g for g in os.listdir() if g.startswith('GEDI02_B') and g.endswith('.h5')] # List all GEDI L2B .h5 files in inDir
gediFiles
Out[3]:
2. Import and Interpret Data ¶
In [4]:
L2B = 'GEDI02_B_2019170155833_O02932_T02267_02_001_01.h5'
L2B
Out[4]:
The standard format for GEDI filenames is as follows:¶
GEDI02_B: 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.¶
In [5]:
gediL2B = h5py.File(L2B, 'r') # Read file using h5py
Navigate the HDF5 file below.¶
In [6]:
list(gediL2B.keys())
Out[6]:
In [7]:
list(gediL2B['METADATA'])
Out[7]:
This contains useful information such as the creation date, PGEVersion, and VersionID. Below, print the file-level metadata attributes.
In [8]:
for g in gediL2B['METADATA']['DatasetIdentification'].attrs: print(g)
In [9]:
print(gediL2B['METADATA']['DatasetIdentification'].attrs['purpose'])
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/.¶
In [10]:
beamNames = [g for g in gediL2B.keys() if g.startswith('BEAM')]
beamNames
Out[10]:
One useful piece of metadata to retrieve from each beam transect is whether it is a full power beam or a coverage beam.¶
In [11]:
for g in gediL2B['BEAM0000'].attrs: print(g)
In [12]:
for b in beamNames:
print(f"{b} is a {gediL2B[b].attrs['description']}")
Below, pick one of the full power beams that will be used to retrieve GEDI L2B shots in Section 3.¶
In [13]:
beamNames = ['BEAM0110']
Identify all the objects in the GEDI HDF5 file below.¶
Note: This step may take a while to complete.
In [14]:
gediL2B_objs = []
gediL2B.visit(gediL2B_objs.append) # Retrieve list of datasets
gediSDS = [o for o in gediL2B_objs if isinstance(gediL2B[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
Out[14]:
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.¶
In [15]:
lonSample, latSample, shotSample, qualitySample, beamSample = [], [], [], [], [] # Set up lists to store data
# Open the SDS
lats = gediL2B[f'{beamNames[0]}/geolocation/lat_lowestmode'][()]
lons = gediL2B[f'{beamNames[0]}/geolocation/lon_lowestmode'][()]
shots = gediL2B[f'{beamNames[0]}/geolocation/shot_number'][()]
quality = gediL2B[f'{beamNames[0]}/l2b_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
Out[15]:
Above is a dataframe containing columns describing the beam, shot number, lat/lon location, and quality information about each shot.¶
In [16]:
# Clean up variables that will no longer be needed
del beamSample, quality, qualitySample, gediL2B_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.¶
In [17]:
# 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.¶
In [18]:
# Convert to a Geodataframe
latslons = gp.GeoDataFrame(latslons)
latslons = latslons.drop(columns=['Latitude','Longitude'])
latslons['geometry']
Out[18]:
Pull out and plot an example shapely point below.¶
In [19]:
latslons['geometry'][0]
Out[19]:
In [20]:
# 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}))
In [21]:
redwoodNP = gp.GeoDataFrame.from_file('RedwoodNP.geojson') # Import geojson as GeoDataFrame
In [22]:
redwoodNP
Out[22]:
In [23]:
redwoodNP['geometry'][0] # Plot GeoDataFrame
Out[23]:
Defining the vdims below will allow you to hover over specific shots and view information about them.¶
In [24]:
# 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
Out[24]:
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.¶
In [25]:
# Call the function for plotting the GEDI points
gv.Polygons(redwoodNP).opts(line_color='red', color=None) * pointVisual(latslons, vdims = vdims)
Out[25]:
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:¶
Side Note: Wondering what the 0's and 1's for l2b_quality_flag mean?¶
In [26]:
print(f"Quality Flag: {gediL2B[b]['l2b_quality_flag'].attrs['description']}")
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 5.1.¶
After finding one of the shots within Redwood NP, find the index for that shot number so that we can find the correct shot 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
In [27]:
del latslons # No longer need the geodataframe used to visualize the full GEDI orbit
4. Work with GEDI L2B Data ¶
The L2B product contains biophysical information derived from the geolocated GEDI return waveforms including total and vertical profiles of canopy cover and Plant Area Index (PAI), the vertical Plant Area Volume Density (PAVD) profile, and Foliage Height Diversity (FHD).¶
Detailed product information can be found on the GEDI L2B Product Page.
4.1 Import and Extract Specific Shots¶
Notice that there are over a thousand datasets available in the GEDI L2B product. In the code blocks below, you will subset to just a few of the datasets available.¶
In this section, learn how to extract and subset specific shots and plot Plant Area Volume Density (PAVD) using holoviews.¶
In [28]:
len(gediSDS)
Out[28]:
In [29]:
beamNames
Out[29]:
In [30]:
beamSDS = [g for g in gediSDS if beamNames[0] in g] # Subset to a single beam
len(beamSDS)
Out[30]:
We will set the shot index used as an example from the GEDI L1B Tutorial and GEDI L2A Tutorial to show how to subset a single shot of GEDI L2B data.¶
In [31]:
shot = 29320619500465599
In [32]:
index = np.where(gediL2B[f'{beamNames[0]}/shot_number'][()]==shot)[0][0] # Set the index for the shot identified above
index
Out[32]:
In [33]:
pavd = gediL2B[[g for g in beamSDS if g.endswith('/pavd_z')][0]] # PAVD
Print the description for the PAVD dataset.¶
In [34]:
print(f"Plant Area Volume Density is {pavd.attrs['description']}")
Below, open the dz layer in order to define the correct vertical step size.¶
In [35]:
# Grab vertical step size
dz = gediL2B[f'{beamNames[0]}/ancillary/dz'][0]
dz
Out[35]:
So the vertical step size is 5.0 meters.¶
In [36]:
print(f"The shape of PAVD is {pavd.shape}.")
In [37]:
# Bring in the desired SDS
elev = gediL2B[f'{beamNames[0]}/geolocation/elev_lowestmode'][()] # Latitude
lats = gediL2B[f'{beamNames[0]}/geolocation/lat_lowestmode'][()] # Latitude
lons = gediL2B[f'{beamNames[0]}/geolocation/lon_lowestmode'][()] # Longitude
Grab the location, elevation, and PAVD metrics for the shot defined above:¶
In [38]:
shotElev = elev[index]
shotLat = lats[index]
shotLon = lons[index]
shotPAVD = pavd[index]
Put everything together to identify the shot that we want to extract:¶
In [39]:
print(f"The shot is located at: {str(shotLat)}, {str(shotLon)} (shot ID: {shot}, index {index}) and is from {beamNames[0]}.")
Next, reformat PAVD into a list of tuples containing each PAVD value and height.¶
In [40]:
pavdAll = []
pavdElev = []
for i, e in enumerate(range(len(shotPAVD))):
if shotPAVD[i] > 0:
pavdElev.append((shot, shotElev + dz * i, shotPAVD[i])) # Append tuple of shot number, elevation, and PAVD
pavdAll.append(pavdElev) # Append to final list
Below, plot each shot by using holoviews Path() function, with the PAVD plotted in the third dimension in shades of green.¶
In [41]:
path1 = hv.Path(pavdAll, vdims='PAVD').options(color='PAVD', clim=(0,0.13), cmap='Greens', line_width=20, colorbar=True,
width=700, height=550, clabel='PAVD', xlabel='Shot Number', tools=['hover'],
ylabel='Elevation (m)', fontsize={'title':16, 'xlabel':16, 'ylabel': 16,
'xticks':12, 'yticks':12,
'clabel':12, 'cticks':10})
path1
Out[41]:
Congratulations! You have plotted your first PAVD profile.¶
In [42]:
# Open all of the desired SDS
dem = gediL2B[[g for g in beamSDS if g.endswith('/digital_elevation_model')][0]][()]
zElevation = gediL2B[[g for g in beamSDS if g.endswith('/elev_lowestmode')][0]][()]
zHigh = gediL2B[[g for g in beamSDS if g.endswith('/elev_highestreturn')][0]][()]
zLat = gediL2B[[g for g in beamSDS if g.endswith('/lat_lowestmode')][0]][()]
zLon = gediL2B[[g for g in beamSDS if g.endswith('/lon_lowestmode')][0]][()]
canopyHeight = gediL2B[[g for g in beamSDS if g.endswith('/rh100')][0]][()]
quality = gediL2B[[g for g in beamSDS if g.endswith('/l2b_quality_flag')][0]][()]
degrade = gediL2B[[g for g in beamSDS if g.endswith('/degrade_flag')][0]][()]
sensitivity = gediL2B[[g for g in beamSDS if g.endswith('/sensitivity')][0]][()]
pavd = gediL2B[f'{beamNames[0]}/pavd_z'][()]
shotNums = gediL2B[f'{beamNames[0]}/shot_number'][()]
# Create a shot index
shotIndex = np.arange(shotNums.size)
In the GEDI L2B product, Canopy Height is stored in units (cm), so below convert to meters.¶
In [43]:
canopyHeight = canopyHeight / 100 # Convert RH100 from cm to m
As mentioned in the sections above, Plant Area Volume Density (pavd) is defined as the Vertical Plant Area Volume Density profile with a vertical step size of dZ. Below, reformat the shape of the PAVD layer in order to add it to the dataframe below.¶
In [44]:
print(f"The shape of Canopy Height is {canopyHeight.shape} vs. the shape of PAVD, which is {pavd.shape}.")
In [45]:
# Set up an empty list to append to
pavdA = []
for i in range(len(pavd)):
# If any of the values are fill value, set to nan
pavdF = [np.nan]
for p in range(len(pavd[i])):
if pavd[i][p]!= -9999:
pavdF.append(pavd[i][p]) # If the value is not fill value, append to list
pavdA.append(pavdF) # Append back to master list
Note: The cell above may take up to a minute to process.¶
Below, notice the reformatted PAVD layer, which should now fit into the dataframe created below.¶
In [46]:
len(pavdA)
Out[46]:
In [47]:
# Take the DEM, GEDI-produced Elevation, and Canopy height 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,
'Plant Area Volume Density': pavdA, 'Sensitivity': sensitivity})
In [48]:
transectDF
Out[48]:
In [49]:
# Plot Canopy Height
canopyVis = hv.Scatter((transectDF['Shot Index'], transectDF['Canopy Height (rh100)']))
canopyVis.opts(color='darkgreen', height=500, width=900, title=f'GEDI L2B Full Transect {beamNames[0]}',
fontsize={'title':16, 'xlabel':16, 'ylabel': 16}, size=0.1, xlabel='Shot Index', ylabel='Canopy Height (m)')
Out[49]: