On-Cloud Data Retrieval and Analysis: Aboveground biomass from GEDI, ICESat-2 and Field Data

Overview¶

This tutorial will demonstrate how to directly access and retrieve the GEDI dataset, ICESat-2 and Field Data in the cloud.

import h5py
import earthaccess
import s3fs
import pandas as pd
import geopandas as gpd
import matplotlib.pyplot as plt
from harmony import BBox, Client, Collection, Request, LinkType
import seaborn as sns
from matplotlib.colors import ListedColormap
sns.set(style='whitegrid')

Area of interest¶

We will use a boundary of the Reserva Florestal Adolpho Ducke, as the region of interest. Reserva Adolpho Ducke is a 10,000-hectare (25,000-acre) protected area of the Amazon rainforest on the outskirts of the city of Manaus, Brazil. It is a part of long term ecological research network and is one of the most intensively studied rainforest in the world.

The boundary file is provided as a GeoJSON file at polygons/reserva_ducke.json. Let’s open and plot this file.

# esri background basemap for maps
xyz = "https://server.arcgisonline.com/ArcGIS/rest/services/World_Imagery/MapServer/tile/{z}/{y}/{x}"
attr = "ESRI"
poly_json = 'polygons/reserva_ducke.json'
poly = gpd.read_file(poly_json) 
poly.explore(color='red', fill=False, tiles=xyz, attr=attr)

Datasets¶

Dubayah, R.O., J. Armston, J.R. Kellner, L. Duncanson, S.P. Healey, P.L. Patterson, S. Hancock, H. Tang, J. Bruening, M.A. Hofton, J.B. Blair, and S.B. Luthcke. 2022. GEDI L4A Footprint Level Aboveground Biomass Density, Version 2.1. ORNL DAAC, Oak Ridge, Tennessee, USA. Dubayah et al. (2022)
Neuenschwander, A. L., Pitts, K. L., Jelley, B. P., Robbins, J., Markel, J., Popescu, S. C., Nelson, R. F., Harding, D., Pederson, D., Klotz, B. & Sheridan, R. (2023). ATLAS/ICESat-2 L3A Land and Vegetation Height. (ATL08, Version 6). Boulder, Colorado USA. NASA National Snow and Ice Data Center Distributed Active Archive Center. Neuenschwander et al. (2023).
dos-Santos, M.N., M.M. Keller, E.R. Pinage, and D.C. Morton. 2022. Forest Inventory and Biophysical Measurements, Brazilian Amazon, 2009-2018. ORNL DAAC, Oak Ridge, Tennessee, USA. Santos et al. (2022)

Workflow¶

Earthdata Authentication¶

We recommend authenticating your Earthdata Login (EDL) information using the earthaccess python library as follows:

auth = earthaccess.login(strategy="netrc") # works if the EDL login already been persisted to a netrc
if not auth.authenticated:
    # ask for EDL credentials and persist them in a .netrc file
    auth = earthaccess.login(strategy="interactive", persist=True)

1. Plot-level Forest Inventory Data¶

1a. Search data files¶

We will use the earthaccess module to search for granules within the dataset Santos et al. (2022).

# CMS Brazil field dataset
doi = '10.3334/ORNLDAAC/2007'
# searching files
granules = earthaccess.search_data(
    doi=doi, # dataset doi
    granule_name = "*.csv",
    bounding_box = tuple(poly.total_bounds.tolist()),
    cloud_hosted = True
)

# print s3 links
[g.data_links(access="direct")[0] for g in granules]

['s3://ornl-cumulus-prod-protected/cms/Forest_Inventory_Brazil/data/DUC_A01_2009_2011_Inventory.csv',
 's3://ornl-cumulus-prod-protected/cms/Forest_Inventory_Brazil/data/DUC_A01_2016_Inventory.csv']

Of the two datasets found, we will retrieve the 2016 inventory which is closer to the GEDI and ICESat-2 temporal range. The Forest Inventory Adolpho Ducke Forest Reserve II (DUC_A01_2016_Inventory) was carried out in Adolpho Ducke Forest Reserve, Amazonas, Brazil. A total of 17 50x50m plots were measured in 2016. Trees with diameter at breast height (DBH) equal to or greater than 35cm were accounted for and measured within the plot area whereas trees with DBH equal to or greater than 10cm were only measured within the subplot area.

1b. Open and plot the data file¶

Let’s open the granule (s3 object) into a xarray. The earthaccess module manages temporary authentication needed for accessing data in NASA’s Earthdata cloud.

for fh in earthaccess.open(granules):
    if fh.granule.data_links()[0].endswith("2016_Inventory.csv"):
        field_df = pd.read_csv(fh, encoding='latin-1')

# printing the dataframe
field_df

# plot tree locations on a map
field_gdf = gpd.GeoDataFrame(field_df, crs="EPSG:4326",
                       geometry=gpd.points_from_xy(field_df["UTM.Northing"],
                                                   field_df["UTM.Easting"]))
m = field_gdf.explore(column='plot', cmap='tab20', legend=None, tiles=xyz, attr=attr)
m

1c. Compute aboveground biomass¶

We will be using the allometry for the tropical forest from Chave et al. 2005 and wood specific gravity for central Amazon obtained from Chave et al. 2006.

# applying allometry
field_df['AGB_kg']=0.0509*0.652*((field_df.DBH)**2)*(field_df.Htotc)
# computing total agb
total_agb = field_df.groupby(['plot'])['AGB_kg'].sum().reset_index()
# computing agbd in Mgha-1
total_agb['AGBD']=total_agb['AGB_kg']/(0.25*1000)
# plotting agbd

fig = plt.figure(figsize=(15,5))
ax1 = fig.add_subplot(1, 3, 1)
sns.histplot(ax=ax1, data=total_agb, y="AGBD", binwidth=20, kde=True, stat='percent')
ax1.set(ylabel="AGBD (Mg/ha)", 
       title=f"Field Plot AGBD estimates")
ax1.set(ylim=(0, 500))
plt.show()

2. GEDI L4A Biomass¶

2a. Define Harmony Request Parameters¶

Let’s create a Harmony Collection object with the concept_id retrieved above. We will also define the GEDI L4A Dubayah et al. (2022) variables of interest and temporal range.

def create_gedi_vars(variables):
    # gedi beams
    beams = ['BEAM0000', 'BEAM0001', 'BEAM0010', 'BEAM0011', 'BEAM0101', 'BEAM0110', 'BEAM1000', 'BEAM1011']
    # combine variables and beam names
    return [f'/{b}/{v}' for b in beams for v in variables]

# gedi variables
variables_l4a = create_gedi_vars(['agbd', 'l4_quality_flag', 'land_cover_data/pft_class'])

# bounding box
b = field_gdf.total_bounds
bounding_box = BBox(w=b[0], s=b[1], e=b[2], n=b[3])

2b. Create and Submit Harmony Request¶

Now, we can create a Harmony request with variables, temporal range, and bounding box and submit the request using the Harmony client object.

doi = '10.3334/ORNLDAAC/2056' # GEDI L4A DOI

# harmony client
harmony_client = Client()
# concept id
concept_id = earthaccess.search_datasets(doi=doi, 
                                         cloud_hosted= True)[0].concept_id()
# define harmony collection
collection = Collection(id=concept_id)
# define harmony request
request = Request(collection=collection, 
              variables=variables_l4a, 
              spatial=bounding_box,
              ignore_errors=True)
# submit harmony request, will return job id
subset_job_id = harmony_client.submit(request)
# retrieve results in-region
results_l4a = harmony_client.result_urls(subset_job_id, 
                                         show_progress=True,
                                         link_type=LinkType.s3)

A temporary S3 Credentials is needed for read-only, same-region (us-west-2), direct access to S3 objects on the Earthdata cloud. We will use the credentials from the harmony_client and pass the credentials to S3Fs class S3FileSystem.

creds = harmony_client.aws_credentials()
s3 = s3fs.S3FileSystem(anon=False,
                       key=creds['aws_access_key_id'],
                       secret=creds['aws_secret_access_key'],
                       token=creds['aws_session_token'])

2c. Read Subset files¶

Let’s direct access the subsetted h5 files and retrieve its values into the pandas dataframe.

def read_gedi_vars(beam):
    """reads through gedi variable hierarchy"""
    col_names = []
    col_val = []
    # read all variables
    for key, value in beam.items():
        # check if the item is a group
        if isinstance(value, h5py.Group):
            # looping through subgroups
            for key2, value2 in value.items():
                col_names.append(key2)
                col_val.append(value2[:].tolist())
        else:
            col_names.append(key)
            col_val.append(value[:].tolist())
    return col_names, col_val


# define an empty pandas dataframe
gedi_df = pd.DataFrame()
# loop through the Harmony results
for s3_url in results_l4a:
    with s3.open(s3_url, mode='rb') as fh:
        with h5py.File(fh) as hf:
            for v in list(hf.keys()):
                if v.startswith('BEAM'):
                    c_n, c_v = read_gedi_vars(hf[v])
                    # Appending to the subset_df dataframe
                    gedi_df = pd.concat([gedi_df, 
                                         pd.DataFrame(map(list, zip(*c_v)), columns=c_n)])

 [ Processing: 100% ] |###################################################| [|]

gedi_df

2d. Quality Filter and Plot¶

We can now quality filter the dataset and only retrieve the good quality shots for trees and shrub cover plant functional types (PFTs).

# remove duplicate columns if any
gedi_df = gedi_df.loc[:,~gedi_df.columns.duplicated()].copy()
# converting to geojson
gedi_gdf = gpd.GeoDataFrame(gedi_df, 
                                geometry=gpd.points_from_xy(gedi_df.lon_lowestmode, 
                                                            gedi_df.lat_lowestmode),
                                crs="EPSG:4326")
# creating mask with good quality shots and trees/shrubs pft class
mask = (gedi_gdf['l4_quality_flag']==1) & (gedi_gdf['pft_class'] <= 5 )
# plotting
gedi_gdf = gedi_gdf[['lat_lowestmode', 'lon_lowestmode', 'pft_class', 'agbd', 'geometry']]
gedi_gdf[mask].explore("agbd", m=m, vmax=500, cmap = "YlGn", alpha=0.5, 
                       radius=10, legend=True)

As we see above, the PFT of the GEDI shots in the area is classed as 2 or “Evergreen Broadleaf Trees”. We will plot the distribution of the AGBD for good quality shots.

fig = plt.figure(figsize=(15,5))
ax1 = fig.add_subplot(1, 3, 1)
sns.histplot(ax=ax1, data=total_agb, y="AGBD", binwidth=20, kde=True, stat='percent')
ax1.set(ylabel="AGBD (Mg / ha)", ylim=(0, 500),
       title="Field Plot AGBD estimates")
ax2 = fig.add_subplot(1, 3, 2)
sns.histplot(ax=ax2, data=gedi_gdf[mask], binwidth=20, y="agbd", kde=True, stat='percent') 
ax2.set(ylabel="AGBD (Mg / ha)", ylim=(0, 500),
       title="GEDI L4A AGBD estimates")
plt.show()

3. ICESat-2 ATL08¶

doi_atl08 = '10.5067/ATLAS/ATL08.006' # ICESat-2 ATL08 DOI

# define harmony collection
collection = Collection(id=earthaccess.search_datasets(doi=doi_atl08,
                                   cloud_hosted= True)[0].concept_id())
# define harmony request
request = Request(collection=collection, 
              spatial=bounding_box,
              ignore_errors=True)

# submit harmony request
subset_job_id = harmony_client.submit(request)
result_atl08 = harmony_client.result_urls(subset_job_id, show_progress=True, 
                                  link_type=LinkType.s3)

Land Segments¶

ATL08 data are grouped into six ground tracks: gt1l, gt1r, gt2l, gt2r, gt3l, and gt3r. The ground tracks with names ending with l are strong beams, and those ending with r are weak beam types. The land_segments group within each ground track contains photon data stored as aggregates of 100 meters.

Within the land_segments group, the geolocation coordinates are provided in latitude and longitude variables. The variable n_seg_ph contains the number of photons in each segment.

Let’s read those variables from the subset files and put them into a geopandas dataframe.

latitude = []
longitude = []
n_seg_ph = []
ph_h = []
classed_pc_flag = []
chm_50 = []
chm_98 = []
for s3_url in result_atl08:
    with s3.open(s3_url, mode='rb') as fh:
        with h5py.File(fh) as hf:
            for var in list(hf.keys()):
                if var.startswith('gt'):
                    latitude.extend(hf[var]['land_segments']['latitude'][:])
                    longitude.extend(hf[var]['land_segments']['longitude'][:])
                    n_seg_ph.extend(hf[var]['land_segments']['n_seg_ph'][:])
                    ph_h.extend(hf[var]['signal_photons']['ph_h'][:])
                    classed_pc_flag.extend(hf[var]['signal_photons']['classed_pc_flag'][:])
                    chm_98.extend(hf[var]['land_segments']['canopy']['h_canopy'][:])
                    chm_v = hf[var]['land_segments']['canopy']['canopy_h_metrics']
                    # chm = chm_v[:, 8]
                    # chm[chm==chm_v.attrs['_FillValue']] = np.nan
                    chm_50.extend(chm_v[:, 8])

 [ Processing:  14% ] |#######                                            | [/]
Job is running with errors.
 [ Processing: 100% ] |###################################################| [|]

# read the lists into a dataframe
atl08_df = pd.DataFrame(list(zip(latitude,longitude,n_seg_ph)), 
                    columns=["latitude", "longitude", "n_seg_ph"])
# converting pandas dataframe into geopandas dataframe
atl08_gdf = gpd.GeoDataFrame(atl08_df, crs="EPSG:4326",
                         geometry=gpd.points_from_xy(atl08_df.longitude, atl08_df.latitude))
# print the first two rows
atl08_gdf.head()

atl08_gdf.explore(m=m)

# pandas dataframe
atl08_df2 = pd.DataFrame(list(zip(chm_50,chm_98)), columns=["chm_50", "chm_98"])
# gedi allometry
atl08_df2["AGBD"]= -134.77015686035156 + (6.653591632843018 * atl08_df2.chm_50) + (6.687118053436279 * atl08_df2.chm_98) 
# predictor_offset
atl08_df2["AGBD"] = atl08_df2["AGBD"]+100
fig = plt.figure(figsize=(15,5))
ax1 = fig.add_subplot(1, 3, 1)
sns.histplot(ax=ax1, data=total_agb, y="AGBD", binwidth=20, kde=True, stat='percent')
ax1.set(ylabel="AGBD (Mg / ha)", ylim=(0, 500),
       title="Field Plot AGBD estimates")
ax2 = fig.add_subplot(1, 3, 2)
sns.histplot(ax=ax2, data=gedi_gdf[mask], binwidth=20, y="agbd", kde=True, stat='percent') 
ax2.set(ylabel="AGBD (Mg / ha)", ylim=(0, 500),
       title="GEDI L4A AGBD estimates")
ax3 = fig.add_subplot(1, 3, 3)
sns.histplot(ax=ax3, data=atl08_df2, binwidth=20, y="AGBD", kde=True, stat='percent') 
ax3.set(ylabel="AGBD (Mg / ha)", ylim=(0, 500),
       title="ICESat-2 ATL08 AGBD estimates")
plt.show()

Photon Classification¶

Photon information is provided in the signal_photons group within each ground track. Let’s plot the photon height (ph_h) and classification (classed_pc_flag) of these four ATL08 files.

# photon classes
classes = ['noise', 'ground', 'canopy', 'top_of_canopy']   
colors = ListedColormap(['gray','blue','green','red'])
plt.figure(figsize=(12, 3))
scatter = plt.scatter(range(len(ph_h)),ph_h,c=classed_pc_flag, cmap=colors, s=1)
plt.legend(handles=scatter.legend_elements()[0], labels=classes)
plt.ylim(-100, 150) 
plt.ylabel("photon height (m)")
plt.show()

References¶

Dubayah, R. O., Armston, J., Kellner, J. R., Duncanson, L., Healey, S. P., Patterson, P. L., Hancock, S., Tang, H., Bruening, J. M., Hofton, M. A., Blair, J. B., & Luthcke, S. B. (2022). GEDI L4A Footprint Level Aboveground Biomass Density, Version 2.1. ORNL Distributed Active Archive Center. 10.3334/ORNLDAAC/2056
Neuenschwander, A., Pitts, K., Jelley, B., Robbins, J., Markel, J., Popescu, S., Nelson, R., Harding, D., Pederson, D., Klotz, B., & Sheridan, R. (2023). ATLAS/ICESat-2 L3A Land and Vegetation Height, Version 6. NASA National Snow. 10.5067/ATLAS/ATL08.006
dos Santos, M. N., Keller, M. M., Pinage, E. R., & Morton, D. C. (2022). Forest Inventory and Biophysical Measurements, Brazilian Amazon, 2009-2018. ORNL Distributed Active Archive Center. 10.3334/ORNLDAAC/2007
Chave, J., Andalo, C., Brown, S., Cairns, M. A., Chambers, J. Q., Eamus, D., Fölster, H., Fromard, F., Higuchi, N., Kira, T., Lescure, J.-P., Nelson, B. W., Ogawa, H., Puig, H., Riéra, B., & Yamakura, T. (2005). Tree allometry and improved estimation of carbon stocks and balance in tropical forests. Oecologia, 145(1), 87–99. 10.1007/s00442-005-0100-x
Chave, J., Muller-Landau, H. C., Baker, T. R., Easdale, T. A., Steege, H. ter, & Webb, C. O. (2006). REGIONAL AND PHYLOGENETIC VARIATION OF WOOD DENSITY ACROSS 2456 NEOTROPICAL TREE SPECIES. Ecological Applications, 16(6), 2356–2367. https://doi.org/10.1890/1051-0761(2006)016[2356:rapvow]2.0.co;2