Getting Started with GEDI L1B Data in Python¶
This tutorial demonstrates how to work with the Geolocated Waveform (GEDI01_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 L1B data to observe GEDI waveforms over Redwood National Park in northern California.
This tutorial will show how to use Python to open GEDI L1B files, visualize the full 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.
Data Used in the Example:¶
- GEDI L1B Geolocated Waveform Data Global Footprint Level - GEDI01_B.001
- 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.
- Science Dataset (SDS) layers:
- /geolocation/latitude_bin0
- /geolocation/longitude_bin0
- /shot_number
- /stale_return_flag
- /geolocation/degrade
- /rx_sample_count
- /rx_sample_start_index
- /rxwaveform
- /geolocation/elevation_bin0
- /geolocation/elevation_lastbin
- /geolocation/digital_elevation_model
- Science Dataset (SDS) layers:
- 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.
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 - Subset and Visualize Waveforms
4.1 Import and Extract Waveforms
4.2 Visualize Waveforms - Quality Filtering
- Plot Profile Transects
6.1 Subset Beam Transects
6.2 Plot Waveform Transects - Spatial Visualization
7.1 Import, Subset, and Quality Filter All Beams
7.2 Spatial Subsetting
7.3 Visualize All Beams: Elevation - 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 L1B observation from June 19, 2019 (orbit 02932). Use the link below to download the file 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 L1B 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, 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.¶
Direct Link to file:¶
- https://e4ftl01.cr.usgs.gov/GEDI/GEDI01_B.001/2019.06.19/GEDI01_B_2019170155833_O02932_T02267_02_003_01.h5
#### Make sure to download the file into theinDirdirectory defined above.
In [3]:
gediFiles = [g for g in os.listdir() if g.startswith('GEDI01_B') and g.endswith('.h5')] # List all GEDI L1B .h5 files in inDir
gediFiles
Out[3]:
2. Import and Interpret Data ¶
In [4]:
L1B = 'GEDI01_B_2019170155833_O02932_T02267_02_003_01.h5'
L1B
Out[4]:
The standard format for GEDI filenames is as follows:¶
GEDI01_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)
003: GOC SDS (software) release number
01: Granule Production Version
Read in a GEDI HDF5 file using the h5py package.¶
In [5]:
gediL1B = h5py.File(L1B, 'r') # Read file using h5py
Navigate the HDF5 file below.¶
In [6]:
list(gediL1B.keys())
Out[6]:
In [7]:
list(gediL1B['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 gediL1B['METADATA']['DatasetIdentification'].attrs: print(g)
In [9]:
print(gediL1B['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 gediL1B.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 gediL1B['BEAM0000'].attrs: print(g)
In [12]:
for b in beamNames:
print(f"{b} is a {gediL1B[b].attrs['description']}")
Below, pick one of the full power beams that will be used to retrieve GEDI L1B waveforms 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]:
gediL1B_objs = []
gediL1B.visit(gediL1B_objs.append) # Retrieve list of datasets
gediSDS = [o for o in gediL1B_objs if isinstance(gediL1B[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, srfSample, degradeSample, beamSample = [], [], [], [], [], [] # Set up lists to store data
# Open the SDS
lats = gediL1B[f'{beamNames[0]}/geolocation/latitude_bin0'][()]
lons = gediL1B[f'{beamNames[0]}/geolocation/longitude_bin0'][()]
shots = gediL1B[f'{beamNames[0]}/shot_number'][()]
srf = gediL1B[f'{beamNames[0]}/stale_return_flag'][()]
degrade = gediL1B[f'{beamNames[0]}/geolocation/degrade'][()]
# 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])
srfSample.append(srf[i])
degradeSample.append(degrade[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,
'Stale Return Flag': srfSample, 'Degrade': degradeSample})
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, degrade, degradeSample, gediL1B_objs, latSample, lats, lonSample, lons, shotSample, shots, srf, srfSample
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 stale_return_flag and degrade mean?¶
In [26]:
print(f"stale_return_flag: {gediL1B[b]['stale_return_flag'].attrs['description']}")
print(f"degrade: {gediL1B[b]['geolocation']['degrade'].attrs['description']}")
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]:
shot = 29320619900465601
In [28]:
index = np.where(gediL1B[f'{beamNames[0]}/shot_number'][()]==shot)[0][0] # Set the index for the shot identified above
index
Out[28]:
In [29]:
del latslons # No longer need the geodataframe used to visualize the full GEDI orbit
4.1 Import and Extract Waveforms¶
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 1 billion 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.¶
In [30]:
# From the SDS list, use list comprehension to find sample_count, sample_start_index, and rxwaveform
sdsCount = gediL1B[[g for g in gediSDS if g.endswith('/rx_sample_count') and beamNames[0] in g][0]]
sdsStart = gediL1B[[g for g in gediSDS if g.endswith('/rx_sample_start_index') and beamNames[0] in g][0]]
sdsWaveform = [g for g in gediSDS if g.endswith('/rxwaveform') and beamNames[0] in g][0]
Print the description for each of these datasets to better understand how we will use them to extract specific waveforms.¶
In [31]:
print(f"rxwaveform is {gediL1B[sdsWaveform].attrs['description']}")
Next, read how to use the rx_sample_count and rx_sample_start_index layers:¶
In [32]:
print(f"rx_sample_count is {sdsCount.attrs['description']}")
print(f"rx_sample_start_index is {sdsStart.attrs['description']}")
Use rx_sample_count and rx_sample_start_index to identify the location of each waveform in rxwaveform.¶
In [33]:
wfCount = sdsCount[index] # Number of samples in the waveform
wfStart = int(sdsStart[index] - 1) # Subtract one because python array indexing starts at 0 not 1
Next grab additional information about the shot, including the unique shot_number, and lat/lon location.¶
In [34]:
wfShot = gediL1B[f'{beamNames[0]}/shot_number'][index]
wfLat = gediL1B[f'{beamNames[0]}/geolocation/latitude_bin0'][index]
wfLon = gediL1B[f'{beamNames[0]}/geolocation/longitude_bin0'][index]
Put everything together to identify the waveform we want to extract:¶
In [35]:
print(f"The waveform located at: {str(wfLat)}, {str(wfLon)} (shot ID: {wfShot}, index {index}) is from beam {beamNames[0]} \
and is stored in rxwaveform beginning at index {wfStart} and ending at index {wfStart + wfCount}")
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.¶
In [36]:
# Grab the elevation recorded at the start and end of the full waveform capture
zStart = gediL1B[f'{beamNames[0]}/geolocation/elevation_bin0'][index] # Height of the start of the rx window
zEnd = gediL1B[f'{beamNames[0]}/geolocation/elevation_lastbin'][index] # Height of the end of the rx window
In [37]:
print("{:,}".format(gediL1B[sdsWaveform].shape[0]))
In [38]:
# Retrieve the waveform sds layer using the sample start index and sample count information to slice the correct dimensions
waveform = gediL1B[sdsWaveform][wfStart: wfStart + wfCount]
In [39]:
# Find elevation difference from start to finish and divide into equal intervals based on sample_count
zStretch = np.add(zEnd, np.multiply(range(wfCount, 0, -1), ((zStart - zEnd) / int(wfCount))))
In [40]:
# match the waveform amplitude values with the elevation and convert to Pandas df
wvDF = pd.DataFrame({'Amplitude (DN)': waveform, 'Elevation (m)': zStretch})
wvDF
Out[40]:
Next, use the Holoviews python package (hv) to plot the waveform.¶
In [41]:
hv.Curve(wvDF) # Basic line graph plotting the waveform
Out[41]:
In [42]:
# Create a holoviews interactive Curve plot with additional parameters defining the plot aesthetics
wfVis = hv.Curve(wvDF).opts(color='darkgreen', tools=['hover'], height=500, width=400,
xlim=(np.min(waveform) - 10, np.max(waveform) + 10), ylim=(np.min(zStretch), np.max(zStretch)),
fontsize={'xticks':10, 'yticks':10,'xlabel':16, 'ylabel': 16, 'title':13}, line_width=2.5, title=f'{str(wfShot)}')
wfVis
Out[42]:
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:¶
Next, plot a couple more waveforms and see if you can capture a shot over the forest.¶
In [43]:
latlons = {} # Set up a dictionary to hold multiple waveforms
latlons[wfShot] = Point(wfLon,wfLat) # Create shapely point and add to dictionary
# Retain waveform and quality layers to be exported later
latlonsWF = [waveform]
latlonsEL = [zStretch]
latlonsSRF = [gediL1B[f'{beamNames[0]}/stale_return_flag'][index]]
latlonsD = [gediL1B[f'{beamNames[0]}/geolocation/degrade'][index]]
In [44]:
# Define which observation to examine by selecting the three indexes prior to the one used above
index = 465597
In [45]:
# Use rx_sample_count and rx_sample_start_index to identify the location of each waveform
wfCount = sdsCount[index] # Number of samples in the waveform
wfStart = int(sdsStart[index] - 1) # Subtract 1, index starts at 0 not 1
wfShot = gediL1B[f'{beamNames[0]}/shot_number'][index] # Unique Shot Number
wfLat = gediL1B[f'{beamNames[0]}/geolocation/latitude_bin0'][index] # Latitude
wfLon = gediL1B[f'{beamNames[0]}/geolocation/longitude_bin0'][index] # Longitude
latlons[wfShot] = Point(wfLon,wfLat) # Create shapely point and add to dictionary
# Grab the elevation recorded at the start and end of the full waveform capture
zStart = gediL1B[f'{beamNames[0]}/geolocation/elevation_bin0'][index] # Height of the start of the rx window
zEnd = gediL1B[f'{beamNames[0]}/geolocation/elevation_lastbin'][index] # Height of the end of the rx window
# Retrieve the waveform sds layer using the sample start index and sample count information to slice the correct dimensions
waveform = gediL1B[sdsWaveform][wfStart: wfStart + wfCount]
In [46]:
# Find elevation difference from start to finish, divide into equal intervals
zStretch = np.add(zEnd, np.multiply(range(wfCount, 0, -1), ((zStart - zEnd) / int(wfCount))))
# match the waveform amplitude values with the elevation and convert to Pandas df
wvDF = pd.DataFrame({'Amplitude (DN)': waveform, 'Elevation (m)': zStretch})
In [47]:
# Append waveform and quality layers to be exported later
latlonsWF.append(waveform)
latlonsEL.append(zStretch)
latlonsSRF.append(gediL1B[f'{beamNames[0]}/stale_return_flag'][index])
latlonsD.append(gediL1B[f'{beamNames[0]}/geolocation/degrade'][index])
In [48]:
# Create a holoviews interactive Curve plot with additional parameters defining the plot aesthetics
visL1B1 = hv.Curve(wvDF).opts(color='darkgreen', tools=['hover'], height=500, width=400,
xlim=(np.min(waveform) - 10, np.max(waveform) + 10), ylim=(np.min(zStretch), np.max(zStretch)),
fontsize={'xticks':10, 'yticks':10,'xlabel':16, 'ylabel': 16, 'title':13}, line_width=2.5, title=f'{str(wfShot)}')
visL1B1
Out[48]:
Now that is starting to look more like a very tall, multi-layered tree canopy! Continue with the next shot:¶
In [49]:
# Define which observation to examine by selecting the three indexes prior to the one used above
index = 465598
In [50]:
# Use rx_sample_count and rx_sample_start_index to identify the location of each waveform
wfCount = sdsCount[index] # Number of samples in the waveform
wfStart = int(sdsStart[index] - 1) # Subtract 1, index starts at 0 not 1
wfShot = gediL1B[f'{beamNames[0]}/shot_number'][index] # Unique Shot Number
wfLat = gediL1B[f'{beamNames[0]}/geolocation/latitude_bin0'][index] # Latitude
wfLon = gediL1B[f'{beamNames[0]}/geolocation/longitude_bin0'][index] # Longitude
latlons[wfShot] = Point(wfLon,wfLat) # Create shapely point and add to dictionary
# Grab the elevation recorded at the start and end of the full waveform capture
zStart = gediL1B[f'{beamNames[0]}/geolocation/elevation_bin0'][index] # Height of the start of the rx window
zEnd = gediL1B[f'{beamNames[0]}/geolocation/elevation_lastbin'][index] # Height of the end of the rx window
# Retrieve the waveform sds layer using the sample start index and sample count information to slice the correct dimensions
waveform = gediL1B[sdsWaveform][wfStart: wfStart + wfCount]
In [51]:
# Find elevation difference from start to finish, divide into equal intervals
zStretch = np.add(zEnd, np.multiply(range(wfCount, 0, -1), ((zStart - zEnd) / int(wfCount))))
# match the waveform amplitude values with the elevation and convert to Pandas df
wvDF = pd.DataFrame({'Amplitude (DN)': waveform, 'Elevation (m)': zStretch})
In [52]:
# Append waveform and quality layers to be exported later
latlonsWF.append(waveform)
latlonsEL.append(zStretch)
latlonsSRF.append(gediL1B[f'{beamNames[0]}/stale_return_flag'][index])
latlonsD.append(gediL1B[f'{beamNames[0]}/geolocation/degrade'][index])
In [53]:
# Create a holoviews interactive Curve plot with additional parameters defining the plot aesthetics
visL1B2 = hv.Curve(wvDF).opts(color='darkgreen', tools=['hover'], height=500, width=400,
xlim=(np.min(waveform) - 10, np.max(waveform) + 10), ylim=(np.min(zStretch), np.max(zStretch)),
fontsize={'xticks':10, 'yticks':10,'xlabel':16, 'ylabel': 16, 'title':13}, line_width=2.5, title=f'{str(wfShot)}')
visL1B2
Out[53]:
In [54]:
# Define which observation to examine by selecting the three indexes prior to the one used above
index = 465599
In [55]:
# Use rx_sample_count and rx_sample_start_index to identify the location of each waveform
wfCount = sdsCount[index] # Number of samples in the waveform
wfStart = int(sdsStart[index] - 1) # Subtract 1, index starts at 0 not 1
wfShot = gediL1B[f'{beamNames[0]}/shot_number'][index] # Unique Shot Number
wfLat = gediL1B[f'{beamNames[0]}/geolocation/latitude_bin0'][index] # Latitude
wfLon = gediL1B[f'{beamNames[0]}/geolocation/longitude_bin0'][index] # Longitude
latlons[wfShot] = Point(wfLon,wfLat) # Create shapely point and add to dictionary
# Grab the elevation recorded at the start and end of the full waveform capture
zStart = gediL1B[f'{beamNames[0]}/geolocation/elevation_bin0'][index] # Height of the start of the rx window
zEnd = gediL1B[f'{beamNames[0]}/geolocation/elevation_lastbin'][index] # Height of the end of the rx window
# Retrieve the waveform sds layer using the sample start index and sample count information to slice the correct dimensions
waveform = gediL1B[sdsWaveform][wfStart: wfStart + wfCount]
In [56]:
# Find elevation difference from start to finish, divide into equal intervals
zStretch = np.add(zEnd, np.multiply(range(wfCount, 0, -1), ((zStart - zEnd) / int(wfCount))))
# match the waveform amplitude values with the elevation and convert to Pandas df
wvDF = pd.DataFrame({'Amplitude (DN)': waveform, 'Elevation (m)': zStretch})
In [57]:
# Append waveform and quality layers to be exported later
latlonsWF.append(waveform)
latlonsEL.append(zStretch)
latlonsSRF.append(gediL1B[f'{beamNames[0]}/stale_return_flag'][index])
latlonsD.append(gediL1B[f'{beamNames[0]}/geolocation/degrade'][index])
In [58]:
# Create a holoviews interactive Curve plot with additional parameters defining the plot aesthetics
visL1B3 = hv.Curve(wvDF).opts(color='darkgreen', tools=['hover'], height=500, width=400,
xlim=(np.min(waveform) - 10, np.max(waveform) + 10), ylim=(np.min(zStretch), np.max(zStretch)),
fontsize={'xticks':10, 'yticks':10,'xlabel':16, 'ylabel': 16, 'title':13}, line_width=2.5, title=f'{str(wfShot)}')
visL1B3
Out[58]:
Lastly, plot the transect of four consecutive waveforms:¶
In [59]:
# The "+" symbol will plot multiple saved holoviews plots together
visL1B1.opts(width=240) + visL1B2.opts(width=240, labelled=[]) + visL1B3.opts(width=240, labelled=[]) + \
wfVis.opts(width=240, labelled=[])
Out[59]:
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.¶
Below, use geoviews to plot the location of the four points to verify what is seen above.¶
In [60]:
# Convert dict to geodataframe
latlons = gp.GeoDataFrame({'Shot Number': list(latlons.keys()),'rxwaveform Amplitude (DN)': latlonsWF,
'rxwaveform Elevation (m)': latlonsEL, 'Stale Return Flag': latlonsSRF, 'Degrade': latlonsD,
'geometry': list(latlons.values())})
In [61]:
latlons['Shot Number'] = latlons['Shot Number'].astype(str) # Convert shot number from integer to string
In [62]:
vdims
Out[62]:
In [63]:
# Create a list of geodataframe columns to be included as attributes in the output map
vdims = []
for f in latlons:
if f not in ['geometry']:
if 'rxwaveform' not in f:
vdims.append(f)
# Plot the geodataframe
gv.Polygons(redwoodNP).opts(line_color='red', color=None) * pointVisual(latlons, vdims=vdims)
Out[63]:
Above, zoom in on the yellow dots until you can get a clearer view of all four shots.¶
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.¶
In [64]:
del wfVis, visL1B1, visL1B2, visL1B3, waveform, wvDF, zStretch, latlonsWF, latlonsEL, latlonsSRF, latlonsD
5. Quality Filtering ¶
Now that you have the desired layers imported as a dataframe for the selected shots, let's perform quality filtering.¶
Below, remove any shots where the stale_return_flag is set to 1 (indicates that a "stale" cue point from the coarse search algorithm is being used) by defining those shots as nan.¶
The syntax of the line below can be read as: in the dataframe, find the rows "where" the stale return 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.¶
In [65]:
latlons = latlons.where(latlons['Stale Return Flag'].ne(1)) # Set any stale returns to NaN
Below, quality filter even further by using the degrade flag (Greater than zero if the shot occurs during a degrade period, zero otherwise).¶
In [66]:
latlons = latlons.where(latlons['Degrade'].ne(1))
Below, drop all of the shots that did not pass the quality filtering standards outlined above from the transectDF.¶
In [67]:
latlons = latlons.dropna() # Drop all of the rows (shots) that did not pass the quality filtering above
In [68]:
print(f"Quality filtering complete, {len(latlons)} high quality shots remaining.")
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.¶
In [69]:
del latlons
6.1 Subset Beam Transects¶
Subset down to a smaller transect centered on the waveforms analyzed in the sections above.¶
In [70]:
print(index)
In [71]:
# Grab 50 points before and after the shot visualized above
start = index - 50
end = index + 50
transectIndex = np.arange(start, end, 1)
6.2 Plot Waveform 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.¶
In [72]:
# Calculate along-track distance
distance = np.arange(0.0, len(transectIndex) * 60, 60) # GEDI Shots are spaced 60 m apart
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.¶
In [73]:
# Create a list of tuples containing Shot number, and each waveform value at height (z)
wfList = []
for s, i in enumerate(transectIndex):
# Basic quality filtering from section 5
if gediL1B['BEAM0110/geolocation/degrade'][i] == 0 and gediL1B['BEAM0110/stale_return_flag'][i] == 0:
zStart = gediL1B['BEAM0110/geolocation/elevation_bin0'][i]
zEnd = gediL1B['BEAM0110/geolocation/elevation_lastbin'][i]
zCount = sdsCount[i]
zStretch = np.add(zEnd, np.multiply(range(zCount, 0, -1), ((zStart - zEnd) / int(zCount))))
waveform = gediL1B[sdsWaveform][sdsStart[i]: sdsStart[i] + zCount]
waves = []
for z, w in enumerate(waveform):
waves.append((distance[s], zStretch[z],w)) # Append Distance (x), waveform elevation (y) and waveform amplitude (z)
wfList.append(waves)
else:
print(f"Shot {s} did not pass quality filtering and will be excluded.")
In [74]:
# Plot the waveform values for the transect with the waveform amplitude being plotted as the 'c' or colormap
l1bVis = hv.Path(wfList, vdims='Amplitude').options(color='Amplitude', clim=(200,405), cmap='Greens', line_width=20, width=950,
height=500, clabel='Amplitude',colorbar=True,
xlabel='Distance Along Transect (m)', ylabel='Elevation (m)',
fontsize={'title':16, 'xlabel':16, 'ylabel': 16, 'xticks':12, 'yticks':12,
'clabel':12, 'cticks':10}, colorbar_position='right',
title='GEDI L1B Waveforms across Redwood National Park: June 19, 2019')
l1bVis
Out[74]: