Python LiDAR – Converting HDF5 LiDAR Data (GEDI) to LAS or LAZ Formats

lidarlidar-full-waveformpythonremote sensing

I downloaded GEDI LiDAR data in hdf5 format. How can I open this file or convert it to the standard LAS format (point cloud)?

When I tried opening the hdf file with QGIS and I got the following error:

Invalid Layer: Raster layer Provider is not valid (provider: gdal, URL:filename)

I also tried opening it with Python:

import pandas as pd
pandas.read_hdf('filename')

I got this error

ValueError: No dataset in HDF5 file

Link to the GEDI datapool.

Best Answer

It's not possible to convert GEDI .h5 file to LAS file as including all data. Because .h5 file includes a lot of information about a point (actually it is a window in GEDI .h5 format, not a point). Also, since LAS file has certain attributes for a point not matching attributes/values in .h5 file, you cannot add all information to LAS file. For example, which value in .h5 file does match Z value in LAS file, elevation_bin0 or elevation_lastbin? Etc.

There is no one way/method to convert any .h5 file to another format. It depends on data structure included in .h5 file. It varies from .h5 to .h5. Therefore you should decide the attributes you would like to use.

Using h5py package, you can easily read h5 file. Firstly, let's examine GEDI .h5 file structure.

import h5py
import pandas as pd
import numpy as np

file_path = "path/to/GEDI01_B_2019108002011_O01959_T03909_02_003_01.h5"
f = h5py.File(file_path, 'r')
print(list(f.keys()))

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

There are 8 BEAMXXXX groups and 1 METADATA group. Now, let's see all datasets in all groups.

def get_h5_structure(f, level=0):
    """    prints structure of hdf5 file    """
    for key in f.keys():
        if isinstance(f[key], h5py._hl.dataset.Dataset):
            print(f"{'  '*level} DATASET: {f[key].name}")
        elif isinstance(f[key], h5py._hl.group.Group):
            print(f"{'  '*level} GROUP: {key, f[key].name}")
            level += 1
            get_h5_structure(f[key], level)
            level -= 1

        if f[key].parent.name == "/":
            print("\n"*2)

get_h5_structure(f)

### OUTPUT: (removed some lines) ###
# GROUP: ('BEAM0000', '/BEAM0000')
#   DATASET: /BEAM0000/all_samples_sum
#   GROUP: ('ancillary', '/BEAM0000/ancillary')
#     DATASET: /BEAM0000/ancillary/master_time_epoch
#     DATASET: /BEAM0000/ancillary/mean_samples
#     DATASET: /BEAM0000/ancillary/smoothing_width
#   DATASET: /BEAM0000/beam
#   DATASET: /BEAM0000/channel
#   DATASET: /BEAM0000/delta_time
#   GROUP: ('geolocation', '/BEAM0000/geolocation')
#     DATASET: /BEAM0000/geolocation/altitude_instrument
#     DATASET: /BEAM0000/geolocation/altitude_instrument_error
#     DATASET: /BEAM0000/geolocation/bounce_time_offset_bin0
#     ...
#
# GROUP: ('BEAM0001', '/BEAM0001')
# ...
#  
# GROUP: ('METADATA', '/METADATA')
#   GROUP: ('DatasetIdentification', '/METADATA/DatasetIdentification')

NOTE: I will use datasets in 'BEAM0000' as an example. For other BEAMXXXX group, you should change group variable.

group = "BEAM0000"

# number_of records
n = f[group]["all_samples_sum"].shape[0]
print(n)
# OUT: 249810

Let's find the keys which have 249810 (n) records. We'll form a DataFrame using these keys. Since there are two nested levels, two for loops are sufficient.

df = pd.DataFrame()

for k, v in f[group].items():
    if isinstance(v, h5py._hl.dataset.Dataset):
        if v.shape[0] == n:
            df[k] = v
    else: # if not dataset, it's group
        # iterate on datasets of the group
        for k2, v2 in v.items():
            if v2.shape[0] == n:
                df[k2] = v2

print(df.head())

### OUTPUT
#    all_samples_sum   beam  channel    delta_time   altitude_instrument  ...  tx_gloc  tx_gloc_error  tx_pulseflag  tx_sample_count  tx_sample_start_index  
# 0         16167838     0        0   4.078333e+07        411250.214378  ...      0.0            0.0             0              128                      1
# 1         16165121     0        0   4.078333e+07        411250.181709  ...      0.0            0.0             0              128                    129
# 2         16180451     0        0   4.078333e+07        411250.149040  ...      0.0            0.0             0              128                    257
# 3         16181775     0        0   4.078333e+07        411250.116372  ...      0.0            0.0             0              128                    385
# 4         16159591     0        0   4.078333e+07        411250.083705  ...      0.0            0.0             0              128                    513
# [5 rows x 77 columns]

surface_type, rxwaveform and txwaveform are missing. As far as I understand, rxwaveform and txwaveform are the most important keys in data.

Let's add surface_type, rxwaveform and txwaveform to df. Please note that each is not a single value, but a list about one point. (See the last 3 columns)

df["surface_type"] = [s_type for s_type in zip(*f[group]["geolocation"]["surface_type"][:,:])]

rxwaveform = f[group]["rxwaveform"][:]
rx_sample_count = f[group]["rx_sample_count"][:]
rx_split_index = f[group]["rx_sample_start_index"][:]-1
df["rxwaveform"] = [ rxwaveform[x:x+i] for x, i in zip(rx_split_index, rx_sample_count) ]

txwaveform = f[group]["txwaveform"][:]
tx_sample_count = f[group]["tx_sample_count"][:]
tx_split_index = f[group]["tx_sample_start_index"][:]-1
df["txwaveform"] = [ txwaveform[x:x+i] for x, i in zip(tx_split_index, tx_sample_count) ]

print(df)

# OUTPUT
#          all_samples_sum  beam  channel    delta_time  altitude_instrument   altitude_instrument_error  .....  tx_pulseflag  tx_sample_count   tx_sample_start_index                  rxwaveform                   txwaveform      surface_type  
#  0              16167838     0        0  4.078333e+07        411250.214378                    0.223205  .....             0              128                       1   [245.59883, 245.52516,...    [246.21742, 246.26566,...   (0, 1, 0, 0, 0)  
#  1              16165121     0        0  4.078333e+07        411250.181709                    0.223205  .....             0              128                     129   [246.6961, 247.62282, ...    [246.30019, 245.81613,...   (0, 1, 0, 0, 0)  
#  ...                 ...   ...      ...           ...                  ...                         ...  .....           ...              ...                     ...                         ...                          ...               ...  
#  249808         16103852     0        0  4.078712e+07        423272.175929                    0.213935  .....             0              128                31975425   [245.15685, 245.5818, ...    [247.31981, 247.07133,...   (0, 1, 0, 0, 0)  
#  249809         16123677     0        0  4.078712e+07        423272.235064                    0.213935  .....             0              128                31975553   [245.64854, 244.94704,...    [247.12299, 247.5319, ...   (0, 1, 0, 0, 0)  
#  
#  [249810 rows x 80 columns]

I don't know what these values mean, therefore, how to use df is up to you.

All necessary script:

import h5py
import pandas as pd
import numpy as np

file_path = "path/to/GEDI01_B_2019108002011_O01959_T03909_02_003_01.h5"
f = h5py.File(file_path, 'r')

group = "BEAM0000"
n = f[group]["all_samples_sum"].shape[0]

df = pd.DataFrame()    
for k, v in f[group].items():
    if isinstance(v, h5py._hl.dataset.Dataset):
        if v.shape[0] == n:
            df[k] = v
    else: # if not dataset, it's group
        # iterate on datasets of the group
        for k2, v2 in v.items():
            if v2.shape[0] == n:
                df[k2] = v2

rxwaveform = f[group]["rxwaveform"][:]
rx_sample_count = f[group]["rx_sample_count"][:]
rx_split_index = f[group]["rx_sample_start_index"][:]-1
df["rxwaveform"] = [ rxwaveform[x:x+i] for x, i in zip(rx_split_index, rx_sample_count)]

txwaveform = f[group]["txwaveform"][:]
tx_sample_count = f[group]["tx_sample_count"][:]
tx_split_index = f[group]["tx_sample_start_index"][:]-1
df["txwaveform"] = [ txwaveform[x:x+i] for x, i in zip(tx_split_index, tx_sample_count)]

df["surface_type"] = [s_type for s_type in zip(*f[group]["geolocation"]["surface_type"][:,:])]

If you prefer, you can save df as shapefile.

import geopandas as gpd

# 2000 sample records
df2 = df[-6000:-4000]

# convert lists to string not to get error
df2['rxwaveform'] = df2['rxwaveform'].apply(str)
df2['txwaveform'] = df2['txwaveform'].apply(str)
df2['surface_type'] = df2['surface_type'].apply(str)

geometries = gpd.points_from_xy(df2.longitude_bin0, df2.latitude_bin0)
gdf = gpd.GeoDataFrame(df2, geometry=geometries)
gdf.crs = '+init=epsg:4326' # WGS84
gdf.to_file("c:/path/to/output.shp")

Best Answer

Related Solutions

Python – Converting Pointcloud Data in UTM from TXT to LAS

Related Question