我最近整理了一点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)
谢谢支持,祝大家工作学习顺利!