Python | 批量裁剪ICESat-2数据(免费代码见文末)

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

for beam in f_in.keys():    if not beam.startswith("gt"):        continue

lat = np.array(in_beam["heights/lat_ph"])lon = np.array(in_beam["heights/lon_ph"])

mask = contains(region, lon, lat)

if np.sum(mask) == 0:    continue

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)

def start(self):    t = threading.Thread(target=self.run)    t.start()

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)

# -*- 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()