【ICESat-2处理】ATL08的读取和可视化（python）

我最近整理了一点ICESat-2的处理代码，我看现在网上的分享内容还挺少的，所以就想把自己整理的发出来吧，方便大家用。我是一个粗人，我不怎么会写代码，不喜欢只分享处理的关键函数，所以我的代码都是从数据读取到最后的可视化都有，直接拿去运行的。请读者们根据自己的需求挑选吧。

前言

其实ATL08的数据信息还挺丰富的，具体的可以去看它的字段介绍，我就照着已有的一些代码，提取出了几种常用的（地面点、冠层点）。如果有需要其它字段，照猫画虎把参数名给改改应该就行了。

因为上一篇ATL03的读取纯纯copy其它博主的，我啥也没改，我就没放原代码，这篇开始就有源代码了，希望能帮到大家。

ATL08数据读取代码

运行前，设置代码中的参数

运行后，我是用窗口选取输入和输出文件，可以根据自己的喜好修改输入输出方式。

【第一步】选择原始的h5文件：

【第二步】设置输出文件

【最后输出一个你设置的研究区内的数据csv表，还有可视化展示】

完整代码

import os
import h5py
import pandas as pd
from tqdm import tqdm
from tkinter import Tk, filedialog
import matplotlib.pyplot as plt


# 设置参数：和ATL03一样，设置三个参数
beams = ['gt3l']  # 待提取的beam的名字'gt1l', 'gt2l', 'gt3l', 'gt1r', 'gt2r',
# 设置经纬度的筛选范围，None 表示不进行筛选
mask_lat = [55.533, 56.550]  # 或 None；练习的话 1°就行了！不然慢
mask_lon = None  # 或 None

# 隐藏主窗口
root = Tk()
root.withdraw()

# 选择输入HDF5文件
file_path = filedialog.askopenfilename(
    title="选择 HDF5 文件",
    filetypes=(("HDF5 files", "*.h5"), ("All files", "*.*"))
)

# 选择输出CSV文件
output_path = filedialog.asksaveasfilename(
    title="选择保存 CSV 文件的位置",
    defaultextension=".csv",
    filetypes=(("CSV files", "*.csv"), ("All files", "*.*"))
)

def plot_data_2(data):
    '''
    可以通过添加多个 ax.scatter 可视化出多个字段信息
    '''
    fig, ax = plt.subplots(figsize=(12, 7))
    ax.scatter(data['latitude'], data['h_canopy_abs'], marker='o', s=5, c='darkgreen', label='ATL08 Abs Canopy')
    ax.scatter(data['latitude'], data['h_surf'], marker='o', s=5, c='sienna', label='ATL08 Ground')
    ax.scatter(data['latitude'], data['h_canopy_rel'], marker='o', s=5, c='green', label='ATL08 Rel Canopy')
    plt.xticks(fontsize=15)
    plt.yticks(fontsize=15)
    ax.set_xlabel('latitude(°N)', fontsize=18)
    ax.set_ylabel('height(m)', fontsize=18)
    plt.rcParams['axes.unicode_minus'] = False
    plt.legend()
    plt.show()



try:
    f = h5py.File(file_path, 'r')  # 打开h5文件
    print('正在处理文件：' + str(file_path))
except Exception as e:
    print(f'读取错误！{e}')
else:
    # 遍历所有需要提取的beam
    for beam in tqdm(beams, desc="处理进度"):
        # 这里因为有些ATL08数据的beam不全，需要跳过
        try:
            gt = f[str(beam)]
            land_segments = gt['land_segments']
        except Exception as e:
            print(f'there is an beam error: {e}')
        else:
            # 获取ATL08数据中需要的属性
            canopy = land_segments['canopy']
            terrain = land_segments['terrain']

            lat = land_segments['latitude'][:]
            lon = land_segments['longitude'][:]

            # 筛选数据，保留经纬度在指定范围内的点
            if mask_lat is not None:
                lat_mask = (lat >= mask_lat[0]) & (lat <= mask_lat[1])
            else:
                lat_mask = [True] * len(lat)

            if mask_lon is not None:
                lon_mask = (lon >= mask_lon[0]) & (lon <= mask_lon[1])
            else:
                lon_mask = [True] * len(lon)

            mask = lat_mask & lon_mask
            if not mask.any():
                continue

            lat = lat[mask]
            lon = lon[mask]
            rgt = land_segments['rgt'][:][mask]
            beam_list = [str(beam) for _ in range(len(lat))]
            cloud_flag_atm = land_segments['cloud_flag_atm'][:][mask]

            time_str = os.path.basename(file_path)[6:14]
            time_nor = [time_str for _ in range(len(lat))]

            h_canopy_rel = canopy['h_canopy'][:][mask]
            h_canopy_abs = canopy['h_canopy_abs'][:][mask]
            h_canopy_uncertainty = canopy['h_canopy_uncertainty'][:][mask]
            h_surf = terrain['h_te_best_fit'][:][mask]

            # 新添加的属性
            h_dif_canopy = canopy['h_dif_canopy'][:][mask]
            h_mean_canopy = canopy['h_mean_canopy'][:][mask]
            h_median_canopy = canopy['h_median_canopy'][:][mask]
            h_min_canopy = canopy['h_min_canopy'][:][mask]
            night_flag = land_segments['night_flag'][:][mask]
            snr = land_segments['snr'][:][mask]

            # 创建新的数据结构来存储数据
            newdata_dic = {
                'latitude': lat,
                'longitude': lon,
                'rgt': rgt,
                'beam': beam_list,
                'cloud_flag': cloud_flag_atm,
                'time': time_nor,
                'h_canopy_rel': h_canopy_rel,
                'h_canopy_abs': h_canopy_abs,
                'h_canopy_uncertainty': h_canopy_uncertainty,
                'h_surf': h_surf,
                'h_dif_canopy': h_dif_canopy,
                'h_mean_canopy': h_mean_canopy,
                'h_median_canopy': h_median_canopy,
                'h_min_canopy': h_min_canopy,
                'night_flag': night_flag,
                'snr': snr,
            }

            newdata = pd.DataFrame(newdata_dic)

            # 设置一些约束：剔除树高大于1000的噪声点
            newdata = newdata[newdata['h_canopy_rel'] <= 1000]

            # 输出csv结果
            newdata.to_csv(output_path, index=False)

            #可视化出来
            plot_data_2(newdata)

谢谢支持，祝大家工作学习顺利！