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

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']

L2B = 'GEDI02_B_2019170155833_O02932_02_T02267_02_003_01_V002.h5'
L2B

'GEDI02_B_2019170155833_O02932_02_T02267_02_003_01_V002.h5'

gediL2B = h5py.File(L2B, 'r')  # Read file using h5py

list(gediL2B.keys())

['BEAM0000',
 'BEAM0001',
 'BEAM0010',
 'BEAM0011',
 'BEAM0101',
 'BEAM0110',
 'BEAM1000',
 'BEAM1011',
 'METADATA']

list(gediL2B['METADATA'])

['DatasetIdentification']

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']

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']

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

# Clean up variables that will no longer be needed
del beamSample, quality, qualitySample, gediL2B_objs, latSample, lats, lonSample, lons, shotSample, shots 

# 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)

# 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

latslons['geometry'][0]

# 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

redwoodNP['geometry'][0]  # Plot GeoDataFrame

# 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

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']

# Call the function for plotting the GEDI points
gv.Polygons(redwoodNP['geometry']).opts(line_color='red', color=None) * pointVisual(latslons, vdims = vdims)

print(f"Quality Flag: {gediL2B[b]['l2b_quality_flag'].attrs['description']}")

Quality Flag: Flag simpilfying selection of most useful data for Level 2B

del latslons  # No longer need the geodataframe used to visualize the full GEDI orbit

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 = 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

# 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).

# 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

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

# 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

# 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

transectDF = transectDF.where(transectDF['Quality Flag'].ne(0))  # Set any poor quality returns to NaN

transectDF

transectDF = transectDF.where(transectDF['Degrade Flag'] < 1)
transectDF = transectDF.where(transectDF['Sensitivity'] > 0.95)

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.

# 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).

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

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']

# 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.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

# 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'])

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)

(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)

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