GEDI_L3 Search and Download

Authors: Anish Bhusal (UAH), Sumant Jha (MSFC/USRA), Jamison French (DevSeed), Aimee Barciauskas (DevSeed), Alex Mandel (DevSeed)

Date: February 8, 2023

Description: In this tutorial, we will search for GEDI L3 data in NASA’s Earthdata Search, and then download a GeoTIFF file from the ORNL DAAC S3. We will then visualize the GeoTIFF file, which consists of canopy heights.

Run This Notebook

To access and run this tutorial within MAAP’s Algorithm Development Environment (ADE), please refer to the “Getting started with the MAAP” section of our documentation.

Disclaimer: it is highly recommended to run a tutorial within MAAP’s ADE, which already includes packages specific to MAAP, such as maap-py. Running the tutorial outside of the MAAP ADE may lead to errors.

About the Data

GEDI L3 Gridded Land Surface Metrics, Version 2

This dataset provides Global Ecosystem Dynamics Investigation (GEDI) Level 3 (L3) gridded mean canopy height, standard deviation of canopy height, mean ground elevation, standard deviation of ground elevation, and counts of laser footprints per 1-km x 1-km grid cells globally within -52 and 52 degrees latitude. GEDI is attached to the International Space Station (ISS) 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.

GEDI L3 data products are gridded by spatially interpolating Level 2 footprint estimates of canopy cover, canopy height, Leaf Area Index (LAI), vertical foliage profile and their uncertainties. Level 2 data contains terrain elevation, canopy height, RH metrics and Leaf Area Index (LAI). The raw waveforms are collected by GEDI system and processed to provide canopy height and profile metrics. These metrics provide easy to use and interpret information about the vertical distribution of canopy material.

Source: GEDI L3 Gridded Land Surface Metrics, Version 2 Data Set Landing Page

Canopy

Source: https://gedi.umd.edu/data/products/

Figure: Sample of GEDI lidar waveform (left). The light brown area under the curve represents the return energy from the canopy, while the dark brown area signifies the return from the underlying typography. The black line is cumulative return energey, starting from the bottom of the ground return (normalized to 0) to the top of canopy (normalized to 1). The diagram on the right shows the distribution of trees that produced the waveform.

Additional Resources

Importing and Installing Packages

This tutorial assumes you’ve all packages installed. If you haven’t already, uncomment the following lines to install these packages.

[ ]:
# !pip install geopandas
# !pip install contextily
# !pip install backoff
# !pip install folium
# !pip install geojsoncontour
[16]:
from maap.maap import MAAP
import pandas as pd
import folium
from rasterio.plot import show
import rasterio
import boto3
import os

After importing necessary packages, the next step is to initialize MAAP constructor using api.maap-project.org as maap_host argument.

[17]:
maap = MAAP(maap_host="api.maap-project.org")

Now, the next step is to seach granules from the CMR. To generate following query, you can use EarthData search feature from MAAP ADE. Refer to this tutorial for more info.

[18]:
results=maap.searchGranule(cmr_host='cmr.earthdata.nasa.gov',concept_id="C2153683336-ORNL_CLOUD", limit=1000)

The above query gives 1000 results by default. The number of necessary results can be changed using limit argument. We can view the GranuleUR from results using:

[19]:
[result['Granule']['GranuleUR'] for result in results]
[19]:
['GEDI_L3_LandSurface_Metrics_V2.GEDI03_elev_lowestmode_stddev_2019108_2020287_002_02.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_counts_2019108_2020287_002_02.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_counts_2019108_2021104_002_02.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_rh100_mean_2019108_2020287_002_02.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_rh100_stddev_2019108_2020287_002_02.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_elev_lowestmode_stddev_2019108_2021104_002_02.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_rh100_stddev_2019108_2021104_002_02.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_elev_lowestmode_mean_2019108_2020287_002_02.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_elev_lowestmode_mean_2019108_2021104_002_02.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_rh100_mean_2019108_2021104_002_02.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_elev_lowestmode_stddev_2019108_2021216_002_02.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_elev_lowestmode_mean_2019108_2021216_002_02.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_rh100_stddev_2019108_2021216_002_02.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_rh100_mean_2019108_2021216_002_02.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_counts_2019108_2021216_002_02.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_elev_lowestmode_mean_2019108_2022019_002_03.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_counts_2019108_2022019_002_03.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_rh100_mean_2019108_2022019_002_03.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_elev_lowestmode_stddev_2019108_2022019_002_03.tif',
 'GEDI_L3_LandSurface_Metrics_V2.GEDI03_rh100_stddev_2019108_2022019_002_03.tif']

Before downloading a particular tif, let’s catch the collection and file name for first item in results:

[20]:
granule_ur=results[0]['Granule']['GranuleUR'].split(".")
collection_name=granule_ur[0]
file_name=granule_ur[1]
[21]:
print(f"collection name: {collection_name} | file_name: {file_name}")
collection name: GEDI_L3_LandSurface_Metrics_V2 | file_name: GEDI03_elev_lowestmode_stddev_2019108_2020287_002_02

Download file from ORNL DAAC S3

To download the file from the source, temporary s3 credentials are required for maap package. You can explicitly request s3_cred_endpoint for the credentials. The code below wraps that request to get credentials and download the file to your workspace.

[22]:
def get_s3_creds(url):
    return maap.aws.earthdata_s3_credentials(url)

def get_s3_client(s3_cred_endpoint):
    creds=get_s3_creds(s3_cred_endpoint)
    boto3_session = boto3.Session(
            aws_access_key_id=creds['accessKeyId'],
            aws_secret_access_key=creds['secretAccessKey'],
            aws_session_token=creds['sessionToken']
    )
    return boto3_session.client("s3")

def download_s3_file(s3, bucket, collection_name, file_name):
    os.makedirs("/projects/gedi_l3", exist_ok=True) # create directories, as necessary
    download_path=f"/projects/gedi_l3/{file_name}.tif"
    s3.download_file(bucket, f"gedi/{collection_name}/data/{file_name}.tif", download_path)
    return download_path
[23]:
s3_cred_endpoint= 'https://data.ornldaac.earthdata.nasa.gov/s3credentials'
s3=get_s3_client(s3_cred_endpoint)
[24]:
bucket="ornl-cumulus-prod-protected"
download_path=download_s3_file(s3, bucket, collection_name, file_name)
download_path
[24]:
'/projects/gedi_l3/GEDI03_elev_lowestmode_stddev_2019108_2020287_002_02.tif'

Now, we have the file in our local workspace. It’s time to visualize it using rasterio package

[Optional] Visualization using Rasterio

The downloaded file is too big to read and visualize directly so we might need to scale it down and view it as a small thumbnail.

[25]:
def show_thumbnail(path):
    src=rasterio.open(path)
    oview = src.overviews(1)[0]
    thumbnail = src.read(1, out_shape=(1, int(src.height // oview), int(src.width // oview)))
    show(thumbnail)
[26]:
show_thumbnail(download_path)
../../_images/science_GEDI_GEDI_L3_27_0.png

[Optional] Overlay Raster Layer on top of Folium Map

To properly visualize the canopy heights, we need to display the TIF image on the map. The TIF image file may be too memory and compute-intensive for the kernel causing the process to exit.

[27]:
# tif=rasterio.open(download_path)
# arr=tif.read()
# bounds=tif.bounds
[28]:
#import numpy as np

# x1,y1,x2,y2=bounds
# bbox=[(bounds.bottom, bounds.left), (bounds.top, bounds.right)]
# m=folium.Map(location=[14.59, 120.98], zoom_start=10)
# img = folium.raster_layers.ImageOverlay(image=np.moveaxis(arr, 0, -1), bounds=bbox, opacity=0.9, interactive=True, cross_origin=False, zindex=1)
# m