Converting GEDI HDF5 LiDAR data to a format readable in GIS

Question

Based on this code here that converts .h5 format into a shapefile, where 'all_samples_sum' feature was replaced with 'digital_elevation_model'. The code I used is based on this code. The error suggests that there missing variables and some features in multidimensions.

The code below only reads features in 1D array, however, I would like it to also include multidimensional arrays and possibly to collect data in all eight beam groups and put them in one dataframe so that I can export as a shapefile.

gediL2B = h5py.File('C:/Username/processed_GEDI02_A_2019144182427_O02530_04_T03532_02_003_01_V002.h5', 'r')

group = "BEAM0000"
no_of_records = f[group]["digital_elevation_model"].shape[0]
print(no_of_records)

#create empty dataframe:
df = pd.DataFrame()
  
for key, value in gediL2B[group].items():
    if isinstance(value, h5py._hl.dataset.Dataset):
        if (len(value.shape) == 1):
            # create a new DataFrame column
            df[key] = value[:]
    else: # if not dataset, it's group
        # iterate on datasets of the group
        for key2, value2 in value.items():
            if isinstance(value2, h5py._hl.dataset.Dataset) and len(value2.shape) == 0:
                # create a new DataFrame column
                df[key2] = value2[:]
print(df.head())

dataset

Error:

KeyError: 'elevs_allmodes_a1'
   
ValueError: Expected a 1D array, got an array with shape (692, 20)

Answer 1

First, I have to say that there are too many datasets that I don't know what it means in the file.

I will try to explain step by step:

• Import necessary packages

  import h5py
  import pandas as pd
  import numpy as np
  import geopandas as gpd
  

• Open h5 file

  file_path = r"C:/path/to/GEDI.h5"
  root = h5py.File(file_path, 'r')
  

• Make an empty DataFrame to collect all data

  all_df = pd.DataFrame()
  

• Define a function to convert all datasets in groups into columns.

  def dataset_to_df(group, n, key_path='', df=None, first_call=False):
      """
      Converts datasets in a group to a DataFrame.
      key_path: used as column name
      """
      # create a dataframe at the first call
      if first_call:
          df = pd.DataFrame()
          # add beam_id column
          df["beam_id"] = pd.Series([group.name[1:]]*n)    
  
      # iterate over the items in current group
      for key, value in group.items():
          col_name =  key if first_call else f"{key_path}/{key}"
  
          # if the item is a group, recursively call this function  
          if isinstance(group[key], h5py.Group):
              dataset_to_df(group[key], n, col_name, df)
  
          # if the key is a dataset, add it to the dataframe as a column
          else:            
              # if the dataset is the same length as the number of beams
              if value.shape[0] == n:                                
                  if value.ndim == 1: # if the dataset is 1D
  
                      # add each dataset as a column to the dataframe               
                      df[col_name] = value
  
                  else: # if the dataset is 2D, add it as each row is list 
                      df[col_name] = [*value[:]]
  
      return df
  

• Skip METADATA

  keys = list(root.keys())[:-1]
  print(keys)
  
  # OUT: 
  # ['BEAM0000', 'BEAM0001', 'BEAM0010', 'BEAM0011',
  #  'BEAM0101', 'BEAM0110', 'BEAM1000', 'BEAM1011']
  

• Iterate over keys and construct the DataFrame containing all BEAMXXXX

  for key in keys:
      print(key)
      n = root[f"{key}/beam"].shape[0] 
  
      # get a dataframe for the current BEAMXXXX group
      df = dataset_to_df(root[key], n, first_call=True)
  
      # merge the current BEAMXXXX dataframe with the previous ones
      all_df = pd.concat([all_df, df], axis=0, ignore_index = True)
  

• Let's check the DataFrame

  print(all_df)
  
  # OUT:
  #
  #        beam_id  beam  channel  degrade_flag    delta_time ...
  # 0     BEAM0000     0        0            70  4.396231e+07 ...
  # 1     BEAM0000     0        0            70  4.396231e+07 ...
  # ...        ...   ...      ...           ...           ... ...
  # 6164  BEAM1011    11        5            70  4.396232e+07 ...
  # 6165  BEAM1011    11        5            70  4.396232e+07 ... 
  
  # [6166 rows x 404 columns]
  

• Column samples: Note the column names used for nested data. Names are like group/dataset_name. 2D values are added as list.

  all_df.iloc[:, [0, 1, 12, 105]]
  
  # OUT:
  #
  #        beam_id     beam   geolocation/elev_highestreturn_a1   geolocation/rh_a1                                
  # 0      BEAM0000    0      832.284973                          [-482, -449, -419, -389, -366, -344, -325, -30...
  # 1      BEAM0000    0      833.631470                          [-460, -426, -396, -366, -336, -314, -292, -26...
  # ...    ...         ...    ...                                 ...                                              
  # 6164   BEAM1011    11     105.038933                          [-906, -801, -722, -651, -584, -520, -460, -40...
  # 6165   BEAM1011    11     41.946384                           [-168, -164, -157, -153, -149, -142, -138, -13...
  # [6166 rows × 4 columns]
  

• Convert the DataFrame to GeoDataFrame: I use geolocation/longitude_1gfit and geolocation/latitude_1gfit for coordinates. Because I don't know what latitude/longitude to use.

  geometries = gpd.points_from_xy(all_df["geolocation/longitude_1gfit"],
                                  all_df["geolocation/latitude_1gfit"])
  
  gdf = gpd.GeoDataFrame(all_df, geometry=geometries)
  gdf.crs = 'EPSG:4326' # WGS84
  

• Plot the GeoDataFrame: (You can see it if you use jupyter notebook)

  gdf.plot(gdf.plot(figsize=(10, 10), markersize=0.1)
  

  

If you use this data by means of geopandas, no problem. But when you want to save it as a geospatial format, a problem arises here. Because 2D data are added as list and when you try to store it in a Shapefile, SQLite, GeoPackage or GeoJSON, the following line throws ValueError: Invalid field type <class 'list'> error. By the way, avoid using shapefiles for GEDI .h5 files. Because, in shapefile, column names are limited to 10 characters.

gdf.to_file("c:/path/to/output.gpkg", layer='LAYER_NAME' driver='GPKG')

One way to solve this problem is to convert the list to string (e.g. [1, 2, 3] -> '[1, 2, 3]'). But when you need to use the data, you have to convert the string back to integer/float list.

To convert list-type column to string, change (in method definition)

df[col_name] = [*value[:]]

to

df[col_name] = pd.Series(value[:].tolist()).apply(lambda x: str(x))

Now you can run:

gdf.to_file("c:/path/to/output.gpkg", layer='LAYER_NAME' driver='GPKG')

---

All code:

import h5py
import pandas as pd
import numpy as np
import geopandas as gpd

file_path = r"C:/path/to/GEDI.h5"
root = h5py.File(file_path, 'r')

all_df = pd.DataFrame()

def dataset_to_df(group, n, key_path='', df=None, first_call=False):
    if first_call:
        df = pd.DataFrame()
        df["beam_id"] = pd.Series([group.name[1:]]*n)    

    for key, value in group.items():
        col_name =  key if first_call else f"{key_path}/{key}"

        if isinstance(group[key], h5py.Group):
            dataset_to_df(group[key], n, col_name, df)
        else:            
            if value.shape[0] == n:                                
                if value.ndim == 1:
                    df[col_name] = value
                else:
                    df[col_name] = pd.Series(value[:].tolist()).apply(lambda x: str(x))

    return df if df else None

keys = list(root.keys())[:-1]

for key in keys:
    n = root[f"{key}/beam"].shape[0] 
    df = dataset_to_df(root[key], n, first_call=True)
    all_df = pd.concat([all_df, df], axis=0, ignore_index = True)
    
geometries = gpd.points_from_xy(all_df["geolocation/longitude_1gfit"],
                                all_df["geolocation/latitude_1gfit"])
gdf = gpd.GeoDataFrame(all_df, geometry=geometries)
gdf.crs = 'EPSG:4326' # WGS84

gdf.to_file("c:/path/to/output.gpkg", layer='LAYER_NAME' driver='GPKG')

In QGIS: