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
+'bokeh', 'matplotlib')
+ gv.extension(import shapely
+import earthaccess
+import warnings
+from shapely.errors import ShapelyDeprecationWarning
+"ignore", category=ShapelyDeprecationWarning) warnings.filterwarnings(
Getting Started with GEDI L2B Version 2 Data in Python
+This tutorial demonstrates how to work with the Canopy Cover and Vertical Profile Metrics (GEDI02_B.002) 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 Version 1 and Version 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 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.
+-
+
- 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 +
+
Data Used in the Example:
+-
+
- GEDI L2B Canopy Cover and Vertical Profile Metrics Data Global Footprint Level - GEDI02_B.002
+
-
+
- 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
+
+ - /selected_l2a_algorithm +
+
+ - 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).
Topics Covered:
+-
+
- Get Started 1.1 Import Packages
+1.2 Set Up the Working Environment and Retrieve Files
+1.3 Authentication
+
+ - Search for GEDI Granules
+
+ - Import and Interpret Data
+3.1 Open a GEDI HDF5 File and Read File Metadata
+3.2 Read SDS Metadata and Subset by Beam
+
+ - Visualize a GEDI Sub-Orbit
+4.1 Subset by Layer and Create a Geodataframe
+4.2 Visualize a Geodataframe
+ - Work with GEDI L2B Data
+5.1 Import and Extract PAVD
+5.2 Visualize PAVD
+
+ - Work with GEDI L2B Beam Transects
+6.1 Quality Filtering
+6.2 Plot Beam Transects
+6.3 Subset Beam Transects
+
+ - Plot Profile Transects
+7.1 Plot PAVD Transects
+
+ - Spatial Visualization
+8.1 Import, Subset, and Quality Filter all Beams
+8.2 Spatial Subsetting
+8.3 Visualize All Beams: Canopy Height, Elevation, and PAI
+
+ - Export Subsets as GeoJSON Files +
Source Code used to Generate this Tutorial:
+The repository containing all of the required files is located at: https://github.com/nasa/GEDI-Data-Resources
+1. Get Started
+1.1 Import Packages
+Import the required packages and set the input/working directory to run this Jupyter Notebook locally.
+1.2 Set Up the Working Environment and Retrieve Files
+The input directory is defined as the current working directory.
+= os.getcwd() # Set input directory to the current working directory
+ inDir
+ os.chdir(inDir)= inDir.rsplit('2023-ssc')[0] + 'shared/2023SSC/'
+ data_dir data_dir
'/home/jovyan/shared/2023SSC/'
+1.3 Authentication
+Login to your NASA Earthdata account and create a .netrc file using the login function from the earthaccess library. If you do not have an Earthdata Account, you can create one here.
+# authenticate
+ earthaccess.login()
EARTHDATA_USERNAME and EARTHDATA_PASSWORD are not set in the current environment, try setting them or use a different strategy (netrc, interactive)
+You're now authenticated with NASA Earthdata Login
+Using token with expiration date: 12/24/2023
+Using .netrc file for EDL
+<earthaccess.auth.Auth at 0x7fa6508d1a20>
+2. Search for GEDI Granules
+Below, earthaccess
library is used to find GEDI L2B V2 Granules for an area of interest and a temporal range. The .h5 file will be downloaded to the data
folder. You will need to download the files in order to execute this tutorial.
+# Define query parameters
+= ["C1908350066-LPDAAC_ECS"]
+ conceptID = (-91.5,18,-89,19) # Lower lon, lower lat, upper lon, upper lat
+ bbx = ('2022-06-01', '2022-06-30')
+ tempRange # search
+= earthaccess.search_data(
+ results # short_name = 'GEDI02_B'
+ = conceptID,
+ concept_id = bbx,
+ bounding_box =tempRange,
+ temporal=500
+ count
+ )# download
+= earthaccess.download(
+ downloaded_files
+ results,= f'{data_dir}', # Update the directory only to avoid overwritting by attendees
+ local_path )
Granules found: 5
+ Getting 5 granules, approx download size: 1.59 GB
+File GEDI02_B_2022154163608_O19677_02_T01836_02_003_01_V002.h5 already downloaded
+File GEDI02_B_2022155015315_O19683_03_T05652_02_003_01_V002.h5 already downloaded
+File GEDI02_B_2022158145957_O19738_02_T03565_02_003_01_V002.h5 already downloaded
+File GEDI02_B_2022162224038_O19805_03_T06723_02_003_02_V002.h5 already downloaded
+File GEDI02_B_2022159001702_O19744_03_T08957_02_003_01_V002.h5 already downloaded
+= [g for g in os.listdir(data_dir) if g.startswith('GEDI02_B') and g.endswith('.h5')] # List all GEDI L2B .h5 files in inDir
+ gediFiles gediFiles
['GEDI02_B_2022155015315_O19683_03_T05652_02_003_01_V002.h5',
+ 'GEDI02_B_2022159001702_O19744_03_T08957_02_003_01_V002.h5',
+ 'GEDI02_B_2022158145957_O19738_02_T03565_02_003_01_V002.h5',
+ 'GEDI02_B_2022154163608_O19677_02_T01836_02_003_01_V002.h5',
+ 'GEDI02_B_2022162224038_O19805_03_T06723_02_003_02_V002.h5']
+3. Import and Interpret Data
+3.1 Open a GEDI HDF5 File and Read File Metadata
+Read the file using h5py
.
+= f'{data_dir}{gediFiles[0]}'
+ L2B L2B
'/home/jovyan/shared/2023SSC/GEDI02_B_2022155015315_O19683_03_T05652_02_003_01_V002.h5'
+The standard format for GEDI Version 2 filenames is as follows:
+++GEDI02_B: Product Short Name
+
+2022159001702: Julian Date and Time of Acquisition (YYYYDDDHHMMSS)
+O19744: Orbit Number
+03: Sub-Orbit Granule Number (1-4)
+T08957: 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
Read in a GEDI HDF5 file using the h5py
package.
+= h5py.File(L2B, 'r') # Read file using h5py gediL2B
The GEDI HDF5 file contains groups in which data and metadata are stored.
+First, the 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, ": ", gediL2B['METADATA']['DatasetIdentification'].attrs[g])
PGEVersion : 003
+VersionID : 01
+abstract : The GEDI L2B standard data product contains precise latitude, longitude, elevation, height, cover and vertical profile metrics for each laser footprint located on the land surface.
+characterSet : utf8
+creationDate : 2022-09-22T17:38:15.690108Z
+credit : The software that generates the L2B product was implemented within the GEDI Science Data Processing System at the NASA Goddard Space Flight Center (GSFC) in Greenbelt, Maryland in collaboration with the Department of Geographical Sciences at the University of Maryland (UMD).
+fileName : GEDI02_B_2022155015315_O19683_03_T05652_02_003_01_V002.h5
+language : eng
+originatorOrganizationName : UMD/GSFC GEDI-SDPS > GEDI Science Data Processing System
+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).
+shortName : GEDI_L2B
+spatialRepresentationType : along-track
+status : onGoing
+topicCategory : geoscientificInformation
+uuid : f7ecc77a-3372-4f53-9131-f110b4bbfb5c
+3.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/.
+= [g for g in gediL2B.keys() if g.startswith('BEAM')]
+ beamNames 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 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
+Below, pick one of the full power beams that will be used to retrieve GEDI L2B shots in next section.
+Identify all the objects in the GEDI HDF5 file below.
+Note: This step may take a while to complete.
+= []
+ gediL2B_objs # Retrieve list of datasets
+ gediL2B.visit(gediL2B_objs.append) = [o for o in gediL2B_objs if isinstance(gediL2B[o], h5py.Dataset)] # Search for relevant SDS inside data file
+ gediSDS # gediSDS
+for i in gediSDS if beamNames[4] in i] # Print the datasets for a selected beam [i
['BEAM0101/algorithmrun_flag',
+ 'BEAM0101/ancillary/dz',
+ 'BEAM0101/ancillary/l2a_alg_count',
+ 'BEAM0101/ancillary/maxheight_cuttoff',
+ 'BEAM0101/ancillary/rg_eg_constraint_center_buffer',
+ 'BEAM0101/ancillary/rg_eg_mpfit_max_func_evals',
+ 'BEAM0101/ancillary/rg_eg_mpfit_maxiters',
+ 'BEAM0101/ancillary/rg_eg_mpfit_tolerance',
+ 'BEAM0101/ancillary/signal_search_buff',
+ 'BEAM0101/ancillary/tx_noise_stddev_multiplier',
+ 'BEAM0101/beam',
+ 'BEAM0101/channel',
+ 'BEAM0101/cover',
+ 'BEAM0101/cover_z',
+ 'BEAM0101/fhd_normal',
+ 'BEAM0101/geolocation/degrade_flag',
+ 'BEAM0101/geolocation/delta_time',
+ 'BEAM0101/geolocation/digital_elevation_model',
+ 'BEAM0101/geolocation/elev_highestreturn',
+ 'BEAM0101/geolocation/elev_lowestmode',
+ 'BEAM0101/geolocation/elevation_bin0',
+ 'BEAM0101/geolocation/elevation_bin0_error',
+ 'BEAM0101/geolocation/elevation_lastbin',
+ 'BEAM0101/geolocation/elevation_lastbin_error',
+ 'BEAM0101/geolocation/height_bin0',
+ 'BEAM0101/geolocation/height_lastbin',
+ 'BEAM0101/geolocation/lat_highestreturn',
+ 'BEAM0101/geolocation/lat_lowestmode',
+ 'BEAM0101/geolocation/latitude_bin0',
+ 'BEAM0101/geolocation/latitude_bin0_error',
+ 'BEAM0101/geolocation/latitude_lastbin',
+ 'BEAM0101/geolocation/latitude_lastbin_error',
+ 'BEAM0101/geolocation/local_beam_azimuth',
+ 'BEAM0101/geolocation/local_beam_elevation',
+ 'BEAM0101/geolocation/lon_highestreturn',
+ 'BEAM0101/geolocation/lon_lowestmode',
+ 'BEAM0101/geolocation/longitude_bin0',
+ 'BEAM0101/geolocation/longitude_bin0_error',
+ 'BEAM0101/geolocation/longitude_lastbin',
+ 'BEAM0101/geolocation/longitude_lastbin_error',
+ 'BEAM0101/geolocation/shot_number',
+ 'BEAM0101/geolocation/solar_azimuth',
+ 'BEAM0101/geolocation/solar_elevation',
+ 'BEAM0101/l2a_quality_flag',
+ 'BEAM0101/l2b_quality_flag',
+ 'BEAM0101/land_cover_data/landsat_treecover',
+ 'BEAM0101/land_cover_data/landsat_water_persistence',
+ 'BEAM0101/land_cover_data/leaf_off_doy',
+ 'BEAM0101/land_cover_data/leaf_off_flag',
+ 'BEAM0101/land_cover_data/leaf_on_cycle',
+ 'BEAM0101/land_cover_data/leaf_on_doy',
+ 'BEAM0101/land_cover_data/modis_nonvegetated',
+ 'BEAM0101/land_cover_data/modis_nonvegetated_sd',
+ 'BEAM0101/land_cover_data/modis_treecover',
+ 'BEAM0101/land_cover_data/modis_treecover_sd',
+ 'BEAM0101/land_cover_data/pft_class',
+ 'BEAM0101/land_cover_data/region_class',
+ 'BEAM0101/land_cover_data/urban_focal_window_size',
+ 'BEAM0101/land_cover_data/urban_proportion',
+ 'BEAM0101/master_frac',
+ 'BEAM0101/master_int',
+ 'BEAM0101/num_detectedmodes',
+ 'BEAM0101/omega',
+ 'BEAM0101/pai',
+ 'BEAM0101/pai_z',
+ 'BEAM0101/pavd_z',
+ 'BEAM0101/pgap_theta',
+ 'BEAM0101/pgap_theta_error',
+ 'BEAM0101/pgap_theta_z',
+ 'BEAM0101/rg',
+ 'BEAM0101/rh100',
+ 'BEAM0101/rhog',
+ 'BEAM0101/rhog_error',
+ 'BEAM0101/rhov',
+ 'BEAM0101/rhov_error',
+ 'BEAM0101/rossg',
+ 'BEAM0101/rv',
+ 'BEAM0101/rx_processing/algorithmrun_flag_a1',
+ 'BEAM0101/rx_processing/algorithmrun_flag_a2',
+ 'BEAM0101/rx_processing/algorithmrun_flag_a3',
+ 'BEAM0101/rx_processing/algorithmrun_flag_a4',
+ 'BEAM0101/rx_processing/algorithmrun_flag_a5',
+ 'BEAM0101/rx_processing/algorithmrun_flag_a6',
+ 'BEAM0101/rx_processing/pgap_theta_a1',
+ 'BEAM0101/rx_processing/pgap_theta_a2',
+ 'BEAM0101/rx_processing/pgap_theta_a3',
+ 'BEAM0101/rx_processing/pgap_theta_a4',
+ 'BEAM0101/rx_processing/pgap_theta_a5',
+ 'BEAM0101/rx_processing/pgap_theta_a6',
+ 'BEAM0101/rx_processing/pgap_theta_error_a1',
+ 'BEAM0101/rx_processing/pgap_theta_error_a2',
+ 'BEAM0101/rx_processing/pgap_theta_error_a3',
+ 'BEAM0101/rx_processing/pgap_theta_error_a4',
+ 'BEAM0101/rx_processing/pgap_theta_error_a5',
+ 'BEAM0101/rx_processing/pgap_theta_error_a6',
+ 'BEAM0101/rx_processing/rg_a1',
+ 'BEAM0101/rx_processing/rg_a2',
+ 'BEAM0101/rx_processing/rg_a3',
+ 'BEAM0101/rx_processing/rg_a4',
+ 'BEAM0101/rx_processing/rg_a5',
+ 'BEAM0101/rx_processing/rg_a6',
+ 'BEAM0101/rx_processing/rg_eg_amplitude_a1',
+ 'BEAM0101/rx_processing/rg_eg_amplitude_a2',
+ 'BEAM0101/rx_processing/rg_eg_amplitude_a3',
+ 'BEAM0101/rx_processing/rg_eg_amplitude_a4',
+ 'BEAM0101/rx_processing/rg_eg_amplitude_a5',
+ 'BEAM0101/rx_processing/rg_eg_amplitude_a6',
+ 'BEAM0101/rx_processing/rg_eg_amplitude_error_a1',
+ 'BEAM0101/rx_processing/rg_eg_amplitude_error_a2',
+ 'BEAM0101/rx_processing/rg_eg_amplitude_error_a3',
+ 'BEAM0101/rx_processing/rg_eg_amplitude_error_a4',
+ 'BEAM0101/rx_processing/rg_eg_amplitude_error_a5',
+ 'BEAM0101/rx_processing/rg_eg_amplitude_error_a6',
+ 'BEAM0101/rx_processing/rg_eg_center_a1',
+ 'BEAM0101/rx_processing/rg_eg_center_a2',
+ 'BEAM0101/rx_processing/rg_eg_center_a3',
+ 'BEAM0101/rx_processing/rg_eg_center_a4',
+ 'BEAM0101/rx_processing/rg_eg_center_a5',
+ 'BEAM0101/rx_processing/rg_eg_center_a6',
+ 'BEAM0101/rx_processing/rg_eg_center_error_a1',
+ 'BEAM0101/rx_processing/rg_eg_center_error_a2',
+ 'BEAM0101/rx_processing/rg_eg_center_error_a3',
+ 'BEAM0101/rx_processing/rg_eg_center_error_a4',
+ 'BEAM0101/rx_processing/rg_eg_center_error_a5',
+ 'BEAM0101/rx_processing/rg_eg_center_error_a6',
+ 'BEAM0101/rx_processing/rg_eg_chisq_a1',
+ 'BEAM0101/rx_processing/rg_eg_chisq_a2',
+ 'BEAM0101/rx_processing/rg_eg_chisq_a3',
+ 'BEAM0101/rx_processing/rg_eg_chisq_a4',
+ 'BEAM0101/rx_processing/rg_eg_chisq_a5',
+ 'BEAM0101/rx_processing/rg_eg_chisq_a6',
+ 'BEAM0101/rx_processing/rg_eg_flag_a1',
+ 'BEAM0101/rx_processing/rg_eg_flag_a2',
+ 'BEAM0101/rx_processing/rg_eg_flag_a3',
+ 'BEAM0101/rx_processing/rg_eg_flag_a4',
+ 'BEAM0101/rx_processing/rg_eg_flag_a5',
+ 'BEAM0101/rx_processing/rg_eg_flag_a6',
+ 'BEAM0101/rx_processing/rg_eg_gamma_a1',
+ 'BEAM0101/rx_processing/rg_eg_gamma_a2',
+ 'BEAM0101/rx_processing/rg_eg_gamma_a3',
+ 'BEAM0101/rx_processing/rg_eg_gamma_a4',
+ 'BEAM0101/rx_processing/rg_eg_gamma_a5',
+ 'BEAM0101/rx_processing/rg_eg_gamma_a6',
+ 'BEAM0101/rx_processing/rg_eg_gamma_error_a1',
+ 'BEAM0101/rx_processing/rg_eg_gamma_error_a2',
+ 'BEAM0101/rx_processing/rg_eg_gamma_error_a3',
+ 'BEAM0101/rx_processing/rg_eg_gamma_error_a4',
+ 'BEAM0101/rx_processing/rg_eg_gamma_error_a5',
+ 'BEAM0101/rx_processing/rg_eg_gamma_error_a6',
+ 'BEAM0101/rx_processing/rg_eg_niter_a1',
+ 'BEAM0101/rx_processing/rg_eg_niter_a2',
+ 'BEAM0101/rx_processing/rg_eg_niter_a3',
+ 'BEAM0101/rx_processing/rg_eg_niter_a4',
+ 'BEAM0101/rx_processing/rg_eg_niter_a5',
+ 'BEAM0101/rx_processing/rg_eg_niter_a6',
+ 'BEAM0101/rx_processing/rg_eg_sigma_a1',
+ 'BEAM0101/rx_processing/rg_eg_sigma_a2',
+ 'BEAM0101/rx_processing/rg_eg_sigma_a3',
+ 'BEAM0101/rx_processing/rg_eg_sigma_a4',
+ 'BEAM0101/rx_processing/rg_eg_sigma_a5',
+ 'BEAM0101/rx_processing/rg_eg_sigma_a6',
+ 'BEAM0101/rx_processing/rg_eg_sigma_error_a1',
+ 'BEAM0101/rx_processing/rg_eg_sigma_error_a2',
+ 'BEAM0101/rx_processing/rg_eg_sigma_error_a3',
+ 'BEAM0101/rx_processing/rg_eg_sigma_error_a4',
+ 'BEAM0101/rx_processing/rg_eg_sigma_error_a5',
+ 'BEAM0101/rx_processing/rg_eg_sigma_error_a6',
+ 'BEAM0101/rx_processing/rg_error_a1',
+ 'BEAM0101/rx_processing/rg_error_a2',
+ 'BEAM0101/rx_processing/rg_error_a3',
+ 'BEAM0101/rx_processing/rg_error_a4',
+ 'BEAM0101/rx_processing/rg_error_a5',
+ 'BEAM0101/rx_processing/rg_error_a6',
+ 'BEAM0101/rx_processing/rv_a1',
+ 'BEAM0101/rx_processing/rv_a2',
+ 'BEAM0101/rx_processing/rv_a3',
+ 'BEAM0101/rx_processing/rv_a4',
+ 'BEAM0101/rx_processing/rv_a5',
+ 'BEAM0101/rx_processing/rv_a6',
+ 'BEAM0101/rx_processing/rx_energy_a1',
+ 'BEAM0101/rx_processing/rx_energy_a2',
+ 'BEAM0101/rx_processing/rx_energy_a3',
+ 'BEAM0101/rx_processing/rx_energy_a4',
+ 'BEAM0101/rx_processing/rx_energy_a5',
+ 'BEAM0101/rx_processing/rx_energy_a6',
+ 'BEAM0101/rx_processing/shot_number',
+ 'BEAM0101/rx_range_highestreturn',
+ 'BEAM0101/rx_sample_count',
+ 'BEAM0101/rx_sample_start_index',
+ 'BEAM0101/selected_l2a_algorithm',
+ 'BEAM0101/selected_mode',
+ 'BEAM0101/selected_mode_flag',
+ 'BEAM0101/selected_rg_algorithm',
+ 'BEAM0101/sensitivity',
+ 'BEAM0101/shot_number',
+ 'BEAM0101/stale_return_flag',
+ 'BEAM0101/surface_flag']
+There are several datasets for each beam. View the GEDI L2B Dictionary for more details. You can also print the description for desired datasets.
+print('pai: ', gediL2B['BEAM0101/pai'].attrs['description'])
+print('pai_z: ', gediL2B['BEAM0101/pai_z'].attrs['description'])
+print('pavd_z: ', gediL2B['BEAM0101/pavd_z'].attrs['description'])
+print('shot_number: ', gediL2B['BEAM0101/geolocation/shot_number'].attrs['description'])
+print('rh100: ', gediL2B['BEAM0101/rh100'].attrs['description'])
+print('Quality Flag: ', gediL2B['BEAM0101/l2b_quality_flag'].attrs['description'])
pai: Total plant area index
+pai_z: Vertical PAI profile from canopy height (z) to ground (z=0) with a vertical step size of dZ, where cover(z > z_max) = 0
+pavd_z: Vertical Plant Area Volume Density profile with a vertical step size of dZ
+shot_number: Unique shot ID.
+rh100: Height above ground of the received waveform signal start (rh[101] from L2A)
+Quality Flag: Flag simpilfying selection of most useful data for Level 2B
+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.
+4. Visualize a GEDI Orbit
+In the section below, import GEDI L2B 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.
+4.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.
+= [], [], [], [], [] # Set up lists to store data
+ lonSample, latSample, shotSample, qualitySample, beamSample
+# Open the SDS
+= gediL2B[f'{beamNames[0]}/geolocation/lat_lowestmode'][()]
+ lats = gediL2B[f'{beamNames[0]}/geolocation/lon_lowestmode'][()]
+ lons = gediL2B[f'{beamNames[0]}/geolocation/shot_number'][()]
+ shots = gediL2B[f'{beamNames[0]}/l2b_quality_flag'][()]
+ quality
+# Take every 100th shot and append to list
+for i in range(len(shots)):
+if i % 100 == 0:
+ str(shots[i]))
+ shotSample.append(
+ lonSample.append(lons[i])
+ latSample.append(lats[i])
+ qualitySample.append(quality[i])0])
+ beamSample.append(beamNames[
+ # Write all of the sample shots to a dataframe
+= pd.DataFrame({'Beam': beamSample, 'Shot Number': shotSample, 'Longitude': lonSample, 'Latitude': latSample,
+ latslons 'Quality Flag': qualitySample})
+ latslons
+ | Beam | +Shot Number | +Longitude | +Latitude | +Quality Flag | +
---|---|---|---|---|---|
0 | +BEAM0000 | +196830000300204328 | +-150.407079 | +51.480311 | +0 | +
1 | +BEAM0000 | +196830000300204428 | +-150.363476 | +51.474715 | +0 | +
2 | +BEAM0000 | +196830000300204528 | +-150.322943 | +51.470677 | +0 | +
3 | +BEAM0000 | +196830000300204628 | +-150.280632 | +51.465678 | +0 | +
4 | +BEAM0000 | +196830000300204728 | +-150.239797 | +51.461442 | +0 | +
... | +... | +... | +... | +... | +... | +
1532 | +BEAM0000 | +196830000300357528 | +-77.574387 | +0.617701 | +0 | +
1533 | +BEAM0000 | +196830000300357628 | +-77.544344 | +0.575149 | +0 | +
1534 | +BEAM0000 | +196830000300357728 | +-77.514584 | +0.532945 | +0 | +
1535 | +BEAM0000 | +196830000300357828 | +-77.484620 | +0.490324 | +0 | +
1536 | +BEAM0000 | +196830000300357928 | +-77.454636 | +0.447691 | +0 | +
1537 rows × 5 columns
+Above is a dataframe containing columns describing the beam, shot number, lat/lon location, and quality information about each shot.
+Side Note: Wondering what the 0’s and 1’s for l2b_quality_flag
mean?
+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.
+# 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.
+# Take the lat/lon dataframe and convert each lat/lon to a shapely point
+'geometry'] = latslons.apply(lambda row: Point(row.Longitude, row.Latitude), axis=1) latslons[
Next, convert to a Geopandas
GeoDataFrame.
+# Convert to a Geodataframe
+= gp.GeoDataFrame(latslons)
+ latslons = latslons.drop(columns=['Latitude','Longitude'])
+ latslons 'geometry'] latslons[
0 POINT (-150.40708 51.48031)
+1 POINT (-150.36348 51.47471)
+2 POINT (-150.32294 51.47068)
+3 POINT (-150.28063 51.46568)
+4 POINT (-150.23980 51.46144)
+ ...
+1532 POINT (-77.57439 0.61770)
+1533 POINT (-77.54434 0.57515)
+1534 POINT (-77.51458 0.53295)
+1535 POINT (-77.48462 0.49032)
+1536 POINT (-77.45464 0.44769)
+Name: geometry, Length: 1537, dtype: geometry
+Pull out and plot an example shapely
point below.
+4.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,
+ ='yellow', fontsize={'xticks': 10, 'yticks': 10,
+ color'xlabel':16, 'ylabel': 16}))
Import a GeoJSON of Reserva de la Biósfera Calakmul
National Park as an additional GeoDataFrame.
+= gp.GeoDataFrame.from_file(f'{data_dir}calakmul.geojson') # Import GeoJSON as GeoDataFrame calakmul
calakmul
+ | geometry | +
---|---|
0 | +POLYGON ((-89.00000 19.00000, -91.50000 19.000... | +
'geometry'][0] # Plot GeoDataFrame calakmul[
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
+'geometry']).opts(line_color='red', color=None) * pointVisual(latslons, vdims = vdims) gv.Polygons(calakmul[
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:
+5. 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.
+Below, only data for one shot is extracted.
+4] beamNames[
'BEAM0101'
+= [g for g in gediSDS if beamNames[4] in g] # Subset to a single beam
+ beamSDS len(beamSDS)
197
+In the GEDI L2B product, Canopy Height is stored in units (cm), so below convert to meters.
+= canopyHeight / 100 # Convert RH100 from cm to m canopyHeight
Below, notice the reformatted PAVD layer, which should now fit into the dataframe created below.
+# Take the DEM, GEDI-produced Elevation, and Canopy height and add to a Pandas dataframe
+= pd.DataFrame({
+ transectDF '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,
+ 'Selected L2A Algorithm': selectedAlgorithmL2A
+ })
ValueError: All arrays must be of the same length
+ transectDF