Published: April 20, 2021
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 Version 1 and Version 2 products. The L1B and L2 GEDI products are archived and distributed in the HDF-EOS5 file format.
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 Version 2 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 Version 2 files, visualize the sub-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.
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.
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
If you prefer to not install Conda, the same setup and dependencies can be achieved by using another package manager such as pip
.
This tutorial uses the GEDI L2B observation from June 19, 2019 (orbit 02932, sub-orbit 02
). Use the links below to download the files directly from the LP DAAC Data Pool:
- https://e4ftl01.cr.usgs.gov/GEDI/GEDI02_B.002/2019.06.19/GEDI02_B_2019170155833_O02932_02_T02267_02_003_01_V002.h5 (0.4 GB)
The repository containing all of the required files is located at: https://git.earthdata.nasa.gov/projects/LPDUR/repos/gedi-v2-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
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')
inDir = os.getcwd() + os.sep # Set input directory to the current working directory
inDir
defined above. You will need to download the file directly from the LP DAAC Data Pool in order to execute this tutorial.inDir
directory defined above.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
['GEDI02_B_2019170155833_O02932_02_T02267_02_003_01_V002.h5']
h5py
.L2B = 'GEDI02_B_2019170155833_O02932_02_T02267_02_003_01_V002.h5'
L2B
'GEDI02_B_2019170155833_O02932_02_T02267_02_003_01_V002.h5'
GEDI02_B: Product Short Name
2019170155833: Julian Date and Time of Acquisition (YYYYDDDHHMMSS)
O02932: Orbit Number
02: Sub-Orbit Granule Number (1-4)
T02267: Track Number (Reference Ground Track)
02: Positioning and Pointing Determination System (PPDS) type (00 is predict, 01 rapid, 02 and higher is final)
003: PGE Version Number
01: Granule Production Version
V002: Product Version
h5py
package.gediL2B = h5py.File(L2B, 'r') # Read file using h5py
list(gediL2B.keys())
['BEAM0000',
'BEAM0001',
'BEAM0010',
'BEAM0011',
'BEAM0101',
'BEAM0110',
'BEAM1000',
'BEAM1011',
'METADATA']
METADATA
group contains the file-level metadata.list(gediL2B['METADATA'])
['DatasetIdentification']
This contains useful information such as the creation date, PGEVersion, and VersionID. Below, print the file-level metadata attributes.
for g in gediL2B['METADATA']['DatasetIdentification'].attrs: print(g)
PGEVersion
VersionID
abstract
characterSet
creationDate
credit
fileName
language
originatorOrganizationName
purpose
shortName
spatialRepresentationType
status
topicCategory
uuid
print(gediL2B['METADATA']['DatasetIdentification'].attrs['purpose'])
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).
beamNames = [g for g in gediL2B.keys() if g.startswith('BEAM')]
beamNames
['BEAM0000',
'BEAM0001',
'BEAM0010',
'BEAM0011',
'BEAM0101',
'BEAM0110',
'BEAM1000',
'BEAM1011']
for g in gediL2B['BEAM0000'].attrs: print(g)
description
wp-l2-l2b_githash
wp-l2-l2b_version
for b in beamNames:
print(f"{b} is a {gediL2B[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
beamNames = ['BEAM0110']
Note: This step may take a while to complete.
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
['BEAM0110/algorithmrun_flag',
'BEAM0110/ancillary/dz',
'BEAM0110/ancillary/l2a_alg_count',
'BEAM0110/ancillary/maxheight_cuttoff',
'BEAM0110/ancillary/rg_eg_constraint_center_buffer',
'BEAM0110/ancillary/rg_eg_mpfit_max_func_evals',
'BEAM0110/ancillary/rg_eg_mpfit_maxiters',
'BEAM0110/ancillary/rg_eg_mpfit_tolerance',
'BEAM0110/ancillary/signal_search_buff',
'BEAM0110/ancillary/tx_noise_stddev_multiplier']
GeoPandas
GeoDataFrame for the beam specified above.lat_lowestmode
and lon_lowestmode
to create a shapely
point for each GEDI shot location.pandas
dataframe.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
Beam | Shot Number | Longitude | Latitude | Quality Flag | |
---|---|---|---|---|---|
0 | BEAM0110 | 29320600200419869 | -142.755692 | 26.923896 | 0 |
1 | BEAM0110 | 29320600200419969 | -142.736567 | 26.943242 | 0 |
2 | BEAM0110 | 29320600200420069 | -142.717433 | 26.962569 | 0 |
3 | BEAM0110 | 29320600200420169 | -142.698294 | 26.981891 | 0 |
4 | BEAM0110 | 29320600200420269 | -142.679136 | 27.001187 | 0 |
... | ... | ... | ... | ... | ... |
1066 | BEAM0110 | 29320600200526469 | -80.198452 | 51.796858 | 0 |
1067 | BEAM0110 | 29320600200526569 | -80.114677 | 51.797029 | 0 |
1068 | BEAM0110 | 29320600200526669 | -80.032173 | 51.797173 | 0 |
1069 | BEAM0110 | 29320600200526769 | -79.948571 | 51.797245 | 0 |
1070 | BEAM0110 | 29320600200526869 | -79.865698 | 51.797246 | 0 |
1071 rows × 5 columns
# Clean up variables that will no longer be needed
del beamSample, quality, qualitySample, gediL2B_objs, latSample, lats, lonSample, lons, shotSample, shots
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)
Geopandas
GeoDataFrame.# Convert to a Geodataframe
latslons = gp.GeoDataFrame(latslons)
latslons = latslons.drop(columns=['Latitude','Longitude'])
latslons['geometry']
0 POINT (-142.75569 26.92390)
1 POINT (-142.73657 26.94324)
2 POINT (-142.71743 26.96257)
3 POINT (-142.69829 26.98189)
4 POINT (-142.67914 27.00119)
...
1066 POINT (-80.19845 51.79686)
1067 POINT (-80.11468 51.79703)
1068 POINT (-80.03217 51.79717)
1069 POINT (-79.94857 51.79725)
1070 POINT (-79.86570 51.79725)
Name: geometry, Length: 1071, dtype: geometry
shapely
point below.latslons['geometry'][0]
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}))
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
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}))
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
# 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']
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)
(HINT: find where the orbit intersects the west coast of the United States)
l2b_quality_flag
mean?print(f"Quality Flag: {gediL2B[b]['l2b_quality_flag'].attrs['description']}")
Quality Flag: Flag simpilfying selection of most useful data for Level 2B
2932: Orbit Number
06: Beam Number
0: Reserved for future use
02: Sub-orbit Granule Number
004: Minor frame number
65601: Shot index
del latslons # No longer need the geodataframe used to visualize the full GEDI orbit
Detailed product information can be found on the GEDI L2B Product Page.
holoviews
.len(gediSDS)
1576
beamNames
['BEAM0110']
beamSDS = [g for g in gediSDS if beamNames[0] in g] # Subset to a single beam
len(beamSDS)
197
shot = 29320600200465599
index = np.where(gediL2B[f'{beamNames[0]}/shot_number'][()]==shot)[0][0] # Set the index for the shot identified above
index
45730
pavd_z
) and begin exploring how to plot them.pavd = gediL2B[[g for g in beamSDS if g.endswith('/pavd_z')][0]] # PAVD
print(f"Plant Area Volume Density is {pavd.attrs['description']}")
Plant Area Volume Density is Vertical Plant Area Volume Density profile with a vertical step size of dZ
dz
layer in order to define the correct vertical step size.# Grab vertical step size
dz = gediL2B[f'{beamNames[0]}/ancillary/dz'][0]
dz
5.0
print(f"The shape of PAVD is {pavd.shape}.")
The shape of PAVD is (107014, 30).
dz
.elev_lowestmode
, lat_lowestmode
and lon_lowestmode
.# 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
shotElev = elev[index]
shotLat = lats[index]
shotLon = lons[index]
shotPAVD = pavd[index]
print(f"The shot is located at: {str(shotLat)}, {str(shotLon)} (shot ID: {shot}, index {index}) and is from {beamNames[0]}.")
The shot is located at: 41.284809604761435, -124.0311657853589 (shot ID: 29320600200465599, index 45730) and is from BEAM0110.
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
holoviews
Path() function, with the PAVD plotted in the third dimension in shades of green.path1 = hv.Path(pavdAll, vdims='PAVD').options(color='PAVD', clim=(0,0.1), cmap='Greens', line_width=20, colorbar=True,
width=700, height=550, clabel='PAVD', xlabel='Shot Number',
ylabel='Elevation (m)', fontsize={'title':16, 'xlabel':16, 'ylabel': 16,
'xticks':12, 'yticks':12,
'clabel':12, 'cticks':10})
path1
pandas
Dataframe to store the arrays.# 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)
canopyHeight = canopyHeight / 100 # Convert RH100 from cm to m
print(f"The shape of Canopy Height is {canopyHeight.shape} vs. the shape of PAVD, which is {pavd.shape}.")
The shape of Canopy Height is (107014,) vs. the shape of PAVD, which is (107014, 30).
# 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
len(pavdA)
107014
# 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})
transectDF
Shot Index | Shot Number | Latitude | Longitude | Tandem-X DEM | Elevation (m) | Canopy Elevation (m) | Canopy Height (rh100) | Quality Flag | Degrade Flag | Plant Area Volume Density | Sensitivity | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | 0 | 29320600200419869 | 26.923896 | -142.755692 | -999999.000000 | 3693.503906 | 3693.503906 | 0.0 | 0 | 0 | [nan] | 5.029931 |
1 | 1 | 29320600200419870 | 26.924089 | -142.755501 | -999999.000000 | 3693.518311 | 3693.518311 | 0.0 | 0 | 0 | [nan] | -4.620663 |
2 | 2 | 29320600200419871 | 26.924283 | -142.755309 | -999999.000000 | 3693.532715 | 3693.532715 | 0.0 | 0 | 0 | [nan] | 6.613796 |
3 | 3 | 29320600200419872 | 26.924477 | -142.755118 | -999999.000000 | 3693.547119 | 3693.547119 | 0.0 | 0 | 0 | [nan] | -0.788731 |
4 | 4 | 29320600200419873 | 26.924671 | -142.754927 | -999999.000000 | 3693.561523 | 3693.561523 | 0.0 | 0 | 0 | [nan] | -64.327515 |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
107009 | 107009 | 29320600200526878 | 51.797246 | -79.858241 | -36.021603 | 4535.905273 | 4535.905273 | 0.0 | 0 | 80 | [nan] | 0.446188 |
107010 | 107010 | 29320600200526879 | 51.797246 | -79.857414 | -999999.000000 | 4543.574219 | 4543.574219 | 0.0 | 0 | 80 | [nan] | 0.709804 |
107011 | 107011 | 29320600200526880 | 51.797245 | -79.856577 | -999999.000000 | 4515.902832 | 4515.902832 | 0.0 | 0 | 80 | [nan] | -0.816940 |
107012 | 107012 | 29320600200526881 | 51.797245 | -79.855748 | -999999.000000 | 4516.102051 | 4516.102051 | 0.0 | 0 | 80 | [nan] | -43.301651 |
107013 | 107013 | 29320600200526882 | 51.797245 | -79.854919 | -41.067085 | 4515.569336 | 4515.569336 | 0.0 | 0 | 80 | [nan] | 0.408280 |
107014 rows × 12 columns
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 L2B Full Transect {beamNames[0]}',
fontsize={'title':16, 'xlabel':16, 'ylabel': 16}, size=0.1, xlabel='Shot Index', ylabel='Canopy Height (m)')
del canopyVis, canopyHeight, degrade, dem, pavd, pavdA, quality, sensitivity, shotIndex, shotNums, zElevation, zHigh, zLat, zLon
l2b_quality_flag
is set to 0 by defining those shots as nan
.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 | Plant Area Volume Density | Sensitivity | |
---|---|---|---|---|---|---|---|---|---|---|---|---|
0 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
1 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
2 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
3 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
4 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
107009 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
107010 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
107011 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
107012 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
107013 | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN | NaN |
107014 rows × 12 columns
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.transectDF = transectDF.where(transectDF['Degrade Flag'] < 1)
transectDF = transectDF.where(transectDF['Sensitivity'] > 0.95)
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, 15364 high quality shots remaining.
holoviews
Scatter(). Combine the Tandem-X derived elevation, the GEDI-derived elevation, and the Canopy Top Elevation in a combined holoviews plot.# 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':14, 'xlabel':16, 'ylabel': 16},
title=f'{beamNames[0]} Full Transect: {L2B.split(".")[0]}')
print(index)
45730
# 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.26960578006781, -124.05874672300072) and ends at (41.29996571902301, -124.00367502601864).
.loc
.transectDF = transectDF.loc[start:end] # Subset the Dataframe to only the selected region of interest over Redwood NP
# 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
pavdAll = []
for j, s in enumerate(transectDF.index):
pavdShot = transectDF['Plant Area Volume Density'][s]
elevShot = transectDF['Elevation (m)'][s]
pavdElev = []
# Remove fill values
if np.isnan(pavdShot).all():
continue
else:
del pavdShot[0]
for i, e in enumerate(range(len(pavdShot))):
if pavdShot[i] > 0:
pavdElev.append((distance[j], elevShot + dz * i, pavdShot[i])) # Append tuple of distance, elevation, and PAVD
pavdAll.append(pavdElev) # Append to final list
canopyElevation = [p[-1][1] for p in pavdAll] # Grab the canopy elevation by selecting the last value in each PAVD
holoviews
Path() function, with the PAVD plotted in the third dimension in shades of green.import warnings
warnings.filterwarnings('ignore')
path1 = hv.Path(pavdAll, vdims='PAVD').options(color='PAVD', clim=(0,0.3), cmap='Greens', line_width=8, colorbar=True,
width=950, height=500, clabel='PAVD', xlabel='Distance Along Transect (m)',
ylabel='Elevation (m)', fontsize={'title':16, 'xlabel':16, 'ylabel': 16,
'xticks':12, 'yticks':12,
'clabel':12, 'cticks':10})
path1
path2 = hv.Curve((distance, transectDF['Elevation (m)']), label='Ground Elevation').options(color='black', line_width=2)
path3 = hv.Curve((distance, canopyElevation), label='Canopy Top Elevation').options(color='grey', line_width=1.5)
# Plot all three together
path = path1 * path2 * path3
path.opts(height=500,width=980, ylim=(min(transectDF['Elevation (m)']) - 5, max(canopyElevation) + 5),
xlabel='Distance Along Transect (m)', ylabel='Elevation (m)', legend_position='bottom_right',
fontsize={'title':15, 'xlabel':15, 'ylabel': 15, 'xticks': 14, 'yticks': 14, 'legend': 14},
title=f'GEDI L2B {beamNames[0]} PAVD over Redwood National Park on June 19, 2019')
del distance, canopyElevation, pavdAll, pavdElev, pavdShot, transectDF
beamNames = [g for g in gediL2B.keys() if g.startswith('BEAM')]
beamNames
['BEAM0000',
'BEAM0001',
'BEAM0010',
'BEAM0011',
'BEAM0101',
'BEAM0110',
'BEAM1000',
'BEAM1011']
pandas
DataFrame.# Set up lists to store data
shotNum, dem, zElevation, zHigh, zLat, zLon, canopyHeight, quality, degrade, sensitivity, pai, beamI = ([] for i in range(12))
# Loop through each beam and open the SDS needed
for b in beamNames:
[shotNum.append(h) for h in gediL2B[[g for g in gediSDS if g.endswith('/shot_number') and b in g][0]][()]]
[dem.append(h) for h in gediL2B[[g for g in gediSDS if g.endswith('/digital_elevation_model') and b in g][0]][()]]
[zElevation.append(h) for h in gediL2B[[g for g in gediSDS if g.endswith('/elev_lowestmode') and b in g][0]][()]]
[zHigh.append(h) for h in gediL2B[[g for g in gediSDS if g.endswith('/elev_highestreturn') and b in g][0]][()]]
[zLat.append(h) for h in gediL2B[[g for g in gediSDS if g.endswith('/lat_lowestmode') and b in g][0]][()]]
[zLon.append(h) for h in gediL2B[[g for g in gediSDS if g.endswith('/lon_lowestmode') and b in g][0]][()]]
[canopyHeight.append(h) for h in gediL2B[[g for g in gediSDS if g.endswith('/rh100') and b in g][0]][()]]
[quality.append(h) for h in gediL2B[[g for g in gediSDS if g.endswith('/l2b_quality_flag') and b in g][0]][()]]
[degrade.append(h) for h in gediL2B[[g for g in gediSDS if g.endswith('/degrade_flag') and b in g][0]][()]]
[sensitivity.append(h) for h in gediL2B[[g for g in gediSDS if g.endswith('/sensitivity') and b in g][0]][()]]
[beamI.append(h) for h in [b] * len(gediL2B[[g for g in gediSDS if g.endswith('/shot_number') and b in g][0]][()])]
[pai.append(h) for h in gediL2B[f'{b}/pai'][()]]
# 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)': canopyHeight,
'Quality Flag': quality, 'Plant Area Index': pai,'Degrade Flag': degrade, 'Sensitivity': sensitivity})
del beamI, canopyHeight, degrade, dem, gediSDS, pai, quality, sensitivity, zElevation, zHigh, zLat, zLon, shotNum
len(allDF)
790135
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
allDF = allDF.where(allDF['Latitude'] > minLat)
allDF = allDF.where(allDF['Latitude'] < maxLat)
allDF = allDF.where(allDF['Longitude'] > minLon)
allDF = allDF.where(allDF['Longitude'] < maxLon)
allDF = allDF.dropna() # Drop shots outside of the ROI
len(allDF)
4477
# Set any poor quality returns to NaN
allDF = allDF.where(allDF['Quality Flag'].ne(0))
allDF = allDF.where(allDF['Degrade Flag'] < 1)
allDF = allDF.where(allDF['Sensitivity'] > 0.95)
allDF = allDF.dropna()
len(allDF)
2979
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'])
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)
allDF['Canopy Height (rh100)'] = allDF['Canopy Height (rh100)'] / 100 # Convert canopy height from cm to 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)
(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)
Plant Area Index
:(gvts.EsriImagery * gv.Points(allDF, vdims=vdims).options(color='Plant Area Index',cmap='Greens', size=3, tools=['hover'],
clim=(0,1), colorbar=True, clabel='m2/m2',
title='GEDI PAI 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)
.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.gediL2B.filename # L2B Filename
'GEDI02_B_2019170155833_O02932_02_T02267_02_003_01_V002.h5'
outName = gediL2B.filename.replace('.h5', '.json') # Create an output file name using the input file name
outName
'GEDI02_B_2019170155833_O02932_02_T02267_02_003_01_V002.json'
allDF.to_file(outName, driver='GeoJSON') # Export to GeoJSON
del allDF