前言：
ICESat-2卫星搭载的ATLAS激光雷达系统，自2018年发射以来已积累海量全球光子数据。其ATL03产品提供了每个光子事件的三维坐标和时间信息，广泛应用于冰川监测、植被高度反演、湖泊水位变化等领域。然而，单轨ATL03数据常覆盖上千公里，研究者通常只关注某个流域或行政区内的小范围数据。传统裁剪方式依赖GIS软件手动操作，面对数十甚至上百个文件时效率极低。
因此，本文分享基于Python的图形化批量裁剪工具，支持单文件和文件夹批处理。
希望各位同学点个关注，点个小赞，这将是更新的动力，不胜感激❥(^_-)
1
数据
本工具处理的数据为ICESat-2卫星的ATL03产品，存储格式为HDF5。每个HDF5文件包含gt1l、gt1r、gt2l、gt2r、gt3l、gt3r六个波束组，每个波束下有land_segments、heights等多个层级，其中heights/lat_ph和heights/lon_ph数据集记录了每个光子事件的纬度和经度，这是空间裁剪的核心依据。裁剪边界由用户提供的SHP文件定义，支持多边形和多部件几何。数据输入支持单个H5文件或包含多个H5文件的文件夹，输出为保持原始HDF5层级结构的裁剪后文件。
2
方法
2.1 坐标参考系统一
裁剪的第一步是读取SHP文件并统一坐标系。ICESat-2数据基于WGS84地理坐标系，如果SHP是投影坐标系，直接计算会导致空间错位。因此，代码使用geopandas读取SHP后，强制转换为EPSG:4326：
gdf = gpd.read_file(shp_path)if gdf.crs is None: raise ValueError("shp没有CRS")gdf = gdf.to_crs("EPSG:4326")
转换后，使用unary_union将SHP中所有几何要素合并为一个整体裁剪区域，确保多部件或不连续的研究区能被统一处理：region = gdf.geometry.unary_union
2.2 遍历波束并提取坐标
接下来打开HDF5文件，遍历每个波束组。由于ICESat-2有强光束和弱光束之分，代码通过判断组名是否以“gt”开头来筛选：
for beam in f_in.keys(): if not beam.startswith("gt"): continue
对每个波束提取光子经纬度坐标。lat_ph和lon_ph是一维数组，每个元素对应一个光子事件的空间位置：
lat = np.array(in_beam["heights/lat_ph"])lon = np.array(in_beam["heights/lon_ph"])
2.3 向量化空间查询生成掩膜
获得坐标数组后，需要判断每个光子是否落在裁剪区域内。这里不使用逐点循环，而是调用shapely的向量化contains函数，一次性对全部坐标完成包含判断，返回布尔掩膜：
mask = contains(region, lon, lat)
该函数底层基于NumPy实现，可高效处理百万级光子点。如果掩膜全为False，说明该波束与SHP无交集，直接跳过：
if np.sum(mask) == 0: continue
2.4 递归复制与智能裁剪
裁剪的核心在于保留HDF5文件的层级结构，同时只对与光子一一对应的数据集应用掩膜。为此，定义了一个递归函数copy_group，遍历原始文件的所有Group和Dataset：
def copy_group(in_grp, out_grp): for key in in_grp.keys(): item = in_grp[key] ifisinstance(item, h5py.Group): sub = out_grp.create_group(key) copy_group(item, sub) elif isinstance(item, h5py.Dataset): data = item[()] if data.shape == lat.shape: data = data[mask] out_grp.create_dataset(key, data=data)
对于Group，递归进入下一层；对于Dataset，判断其形状是否与经纬度数组一致——若一致则应用掩膜裁剪，若不一致则保留原样。这样设计保证了像轨道参数、元数据等非光子级数据不会被错误裁剪。
2.5 图形界面与多线程
为降低使用门槛，工具采用Tkinter构建图形界面。用户通过按钮分别选择H5文件或文件夹、SHP边界和输出路径。为避免大文件处理时界面卡死，裁剪任务被放在子线程中运行：
def start(self): t = threading.Thread(target=self.run) t.start()
日志框实时输出当前处理的波束、坐标范围和裁剪点数，让用户直观掌握进度。
2.6 批量处理逻辑
在运行逻辑中，工具判断用户选择的是单个文件还是文件夹。单文件模式直接调用clip_single函数；批量模式则通过glob搜索所有H5文件，逐个处理并自动按“原名_clipped.h5”命名输出：
if self.h5_path: clip_single(self.h5_path, out, self.shp_path, self.log)elif self.folder_path: files = glob(os.path.join(self.folder_path, "*.h5")) for f in files: out = f.replace(".h5", "_clipped.h5") clip_single(f, out, self.shp_path, self.log)
3
完整代码
看到最后，如果有用请大家点个关注，点个小赞，这将是更新的动力，不胜感激❥(^_-)
# -*- coding: utf-8 -*-"""# -*- coding: utf-8 -*-@Time: 2026/6/5 20:25@Author: LXX@File: 批量裁剪.py@IDE：PyCharm@Motto：ABC(Always Be Coding)"""import tkinter as tkfrom tkinter import filedialog, messageboximport threadingimport osimport h5pyimport numpy as npimport geopandas as gpdfrom shapely.vectorized import containsfrom glob import glob# =========================# 核心裁剪函数# =========================def clip_single(input_h5, output_h5, shp_path, log_func): """ 裁剪单个ATL03 H5文件 """ log_func("读取 shp...") gdf = gpd.read_file(shp_path) if gdf.crs is None: raise ValueError("shp没有CRS") gdf = gdf.to_crs("EPSG:4326") region = gdf.geometry.unary_union log_func(f"shp范围: {region.bounds}") with h5py.File(input_h5, "r") as f_in, h5py.File(output_h5, "w") as f_out: for beam in f_in.keys(): if not beam.startswith("gt"): continue log_func(f"处理 {beam}") in_beam = f_in[beam] out_beam = f_out.create_group(beam) if "heights/lat_ph" not in in_beam: log_func("缺少lat/lon") continue lat = np.array(in_beam["heights/lat_ph"]) lon = np.array(in_beam["heights/lon_ph"]) log_func(f"lat范围 {np.nanmin(lat)} ~ {np.nanmax(lat)}") mask = contains(region, lon, lat) log_func(f"shp内点数: {np.sum(mask)}") if np.sum(mask) == 0: log_func("无交集") continue def copy_group(in_grp, out_grp): """ 递归复制HDF5组，并应用掩码 """ for key in in_grp.keys(): item = in_grp[key] if isinstance(item, h5py.Group): sub = out_grp.create_group(key) copy_group(item, sub) elif isinstance(item, h5py.Dataset): data = item[()] if data.shape == lat.shape: data = data[mask] out_grp.create_dataset(key, data=data) copy_group(in_beam, out_beam) log_func("✔ 完成裁剪")# =========================# GUI 主程序# =========================class App: """ 图形界面应用程序 """ def __init__(self, root): """ 初始化GUI组件 """ self.root = root self.root.title("ICESAT-2 ATL03 裁剪工具") self.h5_path = "" self.folder_path = "" self.shp_path = "" self.out_path = "" # ---- UI ---- tk.Button(root, text="选择H5文件", command=self.select_h5).pack() tk.Button(root, text="选择文件夹(批量)", command=self.select_folder).pack() tk.Button(root, text="选择shp", command=self.select_shp).pack() tk.Button(root, text="选择输出目录", command=self.select_out).pack() tk.Button(root, text="开始裁剪", bg="green", fg="white", command=self.start).pack() self.text = tk.Text(root, height=20) self.text.pack() # ========================= # 选择文件 # ========================= def select_h5(self): """选择单个H5文件""" self.h5_path = filedialog.askopenfilename(filetypes=[("H5","*.h5")]) self.log("H5: " + self.h5_path) def select_folder(self): """选择包含H5文件的文件夹""" self.folder_path = filedialog.askdirectory() self.log("Folder: " + self.folder_path) def select_shp(self): """选择shp文件""" self.shp_path = filedialog.askopenfilename(filetypes=[("shp","*.shp")]) self.log("shp: " + self.shp_path) def select_out(self): """选择输出目录""" self.out_path = filedialog.askdirectory() self.log("Output: " + self.out_path) # ========================= # 日志 # ========================= def log(self, msg): """在文本框中显示日志信息""" self.text.insert(tk.END, msg + "\n") self.text.see(tk.END) # ========================= # 启动线程 # ========================= def start(self): """启动处理线程""" t = threading.Thread(target=self.run) t.start() # ========================= # 运行逻辑 # ========================= def run(self): """ 执行裁剪的主逻辑 """ if not self.shp_path: messagebox.showerror("错误", "请选择shp") return # ---- 单文件 ---- if self.h5_path: out = self.h5_path.replace(".h5", "_clipped.h5") self.log("开始单文件裁剪...") clip_single(self.h5_path, out, self.shp_path, self.log) # ---- 批量 ---- elif self.folder_path: files = glob(os.path.join(self.folder_path, "*.h5")) self.log(f"批量文件: {len(files)}") for f in files: out = f.replace(".h5", "_clipped.h5") self.log("处理: " + f) clip_single(f, out, self.shp_path, self.log) else: messagebox.showerror("错误", "请选择H5或文件夹")# =========================# 启动# =========================if __name__ == "__main__": root = tk.Tk() app = App(root) root.mainloop()
10个月宝宝每天需要喝多少奶粉？