你提出的需求本质上是：**对一批 TIFF 格式的小时降水栅格数据（同一目录下），使用普通克里金（Ordinary Kriging）插值重采样至两个目标地理分辨率（0.1° 和 0.05°），插值模型采用高斯（Gaussian）半变异函数模型，其余参数使用默认设置**。注意：你原流程中“栅格→点→克里金→栅格”是 ArcMap 的典型 workaround（因 ArcGIS 不支持直接对栅格进行空间插值重采样），但**在 Python 中，我们完全可以跳过“转点”这一步——因为克里金本质是对空间点数据建模；而你的原始 TIFF 是规则网格降水场，它本身不是观测点数据，因此必须先提取有效像元为采样点（即：将非空值的行列坐标 + 值 → (x, y, z) 点集），再以此点集驱动 Ordinary Kriging 插值到新网格上**。 ✅ 关键澄清： - 普通克里金（Ordinary Kriging）**必须输入离散观测点**（x, y, value），不能直接对输入栅格做“克里金重采样”。所以“TIFF → 有效点集 → OK 插值到新网格”是正确且必要的流程。 - “掩膜提取”在你场景中可省略（因范围一致），但为鲁棒性，我们仍建议用 `rasterio` 读取时自动识别 nodata，并仅提取有效值点。 - 分辨率单位为度（°），需确保输入 TIFF 的 CRS 是地理坐标系（如 WGS84, EPSG:4326），否则 0.1° 无意义。我们将显式检查并强制重投影（若需要）。 - 高斯模型（`variogram_model='gaussian'`）由 `scikit-gstat` 或 `pykrige` 支持；这里推荐 **`pykrige`**（专为地理克里金设计，API 清晰，支持 2D/3D，与 `rasterio`/`numpy` 生态无缝集成）。 📌 技术栈选择： - `rasterio`: 安全读写 GeoTIFF，获取 transform、crs、nodata、bounds - `pykrige`: 提供 `OrdinaryKriging`，支持 `gaussian` variogram，可批量预测网格 - `numpy` / `xarray`: 数值处理与结果组织 - `tqdm`: 进度条（大数据友好） - `os` / `glob`: 批量遍历 TIFF 文件 --- ### ✅ 完整可运行 `.ipynb` 兼容代码（含详细注释） > 💡 使用前请 `pip install rasterio pykrige scikit-gstat numpy tqdm matplotlib`（`scikit-gstat` 是 `pykrige` 可选依赖，但推荐安装以支持更多变差函数诊断） ```python # %% [markdown] # # 🌧️ 批量小时降水 TIFF 普通克里金重采样（0.1° & 0.05°） # - 输入：同一目录下所有 `.tif` 降水栅格（WGS84, 单波段） # - 步骤：读取 → 提取有效点 (x,y,z) → Ordinary Kriging (Gaussian) → 插值到目标分辨率网格 → 输出新 TIFF # - 输出：每个输入文件生成两个输出：`{name}_res0.1.tif`, `{name}_res0.05.tif` # - 注意：自动检测 CRS，非 WGS84 将重投影（警告提示） # %% [code] import os import glob import numpy as np import rasterio from rasterio.transform import from_origin from rasterio.crs import CRS from pykrige.ok import OrdinaryKriging from pykrige.uk import UniversalKriging from tqdm import tqdm import warnings warnings.filterwarnings("ignore", category=FutureWarning) # pykrige 旧版警告 # ------------------------------- # 🔧 用户配置区（请按需修改） # ------------------------------- INPUT_DIR = r"./precip_tiff" # ✅ 替换为你存放 .tif 的绝对或相对路径 OUTPUT_DIR = r"./kriged_output" # ✅ 输出目录（会自动创建） TARGET_RESOLUTIONS = [0.1, 0.05] # ✅ 目标分辨率（单位：度，WGS84） VARIOMODEL = "gaussian" # ✅ 普通克里金变差函数模型（可选："linear", "power", "gaussian", "spherical", "exponential"） N_MAX_POINTS = 10000 # ⚠️ 安全限制：单文件最多采样点数（防内存溢出，可调大） # ------------------------------- # 创建输出目录 os.makedirs(OUTPUT_DIR, exist_ok=True) # 获取所有 TIFF 文件 tiff_files = sorted(glob.glob(os.path.join(INPUT_DIR, "*.tif")) + glob.glob(os.path.join(INPUT_DIR, "*.TIF"))) if not tiff_files: raise FileNotFoundError(f"❌ 在 {INPUT_DIR} 中未找到任何 .tif 文件") print(f"✅ 发现 {len(tiff_files)} 个 TIFF 文件，开始处理...\n") # %% [code] def extract_valid_points(src_ds): """ 从 rasterio dataset 中提取所有有效值点 (lon, lat, value) 自动处理 CRS 转换：确保输出为 WGS84 经纬度 """ # 读取数据和地理信息 data = src_ds.read(1).astype(np.float64) # 强制 float64 避免 pykrige 类型问题 nodata = src_ds.nodata if nodata is not None: data = np.where(data == nodata, np.nan, data) # 获取地理变换（affine transform）和 CRS transform = src_ds.transform crs = src_ds.crs # 计算所有像素中心坐标（行、列 → 地理坐标） rows, cols = np.where(np.isfinite(data)) # 仅取有限值（非 NaN, 非 nodata） if len(rows) == 0: raise ValueError("❌ 栅格中无有效值（全为 nodata 或 NaN）") # 向量化计算地理坐标：transform * (col+0.5, row+0.5) xs, ys = rasterio.transform.xy(transform, rows, cols, offset='center') xs, ys = np.array(xs), np.array(ys) # 若 CRS 不是 WGS84，则重投影到 EPSG:4326 if crs is None or crs.to_epsg() != 4326: from rasterio.warp import reproject, calculate_default_transform, Resampling from rasterio.crs import CRS # 构造临时内存 dataset 用于重投影坐标（更高效做法是用 pyproj） try: import pyproj transformer = pyproj.Transformer.from_crs( crs, CRS.from_epsg(4326), always_xy=True ) lons, lats = transformer.transform(xs, ys) except Exception as e: raise RuntimeError(f"❌ 坐标转换失败，请确认输入栅格有有效 CRS：{e}") else: lons, lats = xs, ys values = data[rows, cols] # 采样点数量控制（随机下采样，保持空间代表性） if len(lons) > N_MAX_POINTS: indices = np.random.choice(len(lons), N_MAX_POINTS, replace=False) lons, lats, values = lons[indices], lats[indices], values[indices] print(f" ⚠️ 点数超限({len(lons)}→{N_MAX_POINTS})，已随机下采样") return lons, lats, values def create_target_grid(bounds, res_deg): """根据输入栅格 bounds 和目标分辨率，生成规则经纬度网格""" left, bottom, right, top = bounds # 向外扩展半像元，确保覆盖完整范围（避免边界截断） left -= res_deg * 0.5 right += res_deg * 0.5 bottom -= res_deg * 0.5 top += res_deg * 0.5 # 生成一维经纬度向量 lon_vec = np.arange(left, right + res_deg * 0.1, res_deg) # +0.1 防止浮点误差漏末尾 lat_vec = np.arange(bottom, top + res_deg * 0.1, res_deg)[::-1] # lat 从北到南（rasterio 标准） # 生成网格（lat 优先，符合 rasterio 行优先） lon_grid, lat_grid = np.meshgrid(lon_vec, lat_vec) return lon_grid, lat_grid, lon_vec, lat_vec def kriging_to_grid(lons, lats, values, lon_grid, lat_grid, variomodel="gaussian"): """ 使用 OrdinaryKriging 对网格插值 返回：2D 插值结果数组（shape = lat_grid.shape） """ # 初始化 OK（自动拟合变差函数） OK = OrdinaryKriging( lons, lats, values, variogram_model=variomodel, verbose=False, enable_plotting=False, nlags=10, # 默认分组数 weight=True # 启用距离加权 ) # 执行批量预测（flatten 后再 reshape） pred_flat = OK.execute("grid", lon_grid[0, :], lat_grid[:, 0]) # 注意：pykrige 的 grid 返回 (y_size, x_size)，与 lon_grid/lat_grid 一致 return pred_flat def write_geotiff(filename, data, transform, crs, nodata=np.nan): """将 2D 数组写为带地理参考的 GeoTIFF""" with rasterio.open( filename, 'w', driver='GTiff', height=data.shape[0], width=data.shape[1], count=1, dtype=data.dtype, crs=crs, transform=transform, nodata=nodata, compress='lzw' ) as dst: dst.write(data, 1) # %% [code] # 主循环：逐文件处理 for tif_path in tqdm(tiff_files, desc="📦 处理文件"): try: base_name = os.path.splitext(os.path.basename(tif_path))[0] print(f"\n🔍 正在处理：{base_name}") # 1️⃣ 读取输入栅格 with rasterio.open(tif_path) as src: # 检查 CRS if src.crs is None: print(f" ⚠️ 输入栅格无 CRS，假设为 EPSG:4326（WGS84）") crs_out = CRS.from_epsg(4326) else: crs_out = src.crs # 2️⃣ 提取有效点 (lon, lat, precip) lons, lats, values = extract_valid_points(src) print(f" ✅ 提取 {len(lons)} 个有效点") # 3️⃣ 为每个目标分辨率构建输出网格 for res in TARGET_RESOLUTIONS: print(f" 🎯 开始插值至 {res}° 分辨率...") # 创建目标网格（经纬度） lon_grid, lat_grid, lon_vec, lat_vec = create_target_grid(src.bounds, res) # 4️⃣ 普通克里金插值 pred_data = kriging_to_grid(lons, lats, values, lon_grid, lat_grid, VARIOMODEL) # 5️⃣ 构建输出 transform（左上角 + 像元大小） # transform = from_origin(west, north, xsize, ysize) transform_out = from_origin( lon_vec[0], # west = min lon lat_vec[0], # north = max lat（因 lat_vec 已[::-1]，首元素是北界） res, res ) # 6️⃣ 写入 GeoTIFF out_name = f"{base_name}_res{res:.2f}.tif".replace(".", "p") # 避免文件名小数点问题 out_path = os.path.join(OUTPUT_DIR, out_name) write_geotiff(out_path, pred_data, transform_out, crs_out, nodata=np.nan) print(f" ✅ 已保存：{out_name}") print(f"✅ 完成 {base_name}") except Exception as e: print(f"❌ 处理 {os.path.basename(tif_path)} 时出错：{e}") continue print(f"\n🎉 全部完成！结果保存在：{OUTPUT_DIR}") ``` --- ### 🔍 代码关键解释： | 模块 | 说明 | |------|------| | **`extract_valid_points()`** | 核心预处理：读取 TIFF → 获取所有非 nodata 像元的地理坐标（自动重投影到 WGS84）→ 返回 `(lons, lats, values)`。若点过多（>10000），随机下采样保证效率。 | | **`create_target_grid()`** | 根据输入栅格 `bounds`（左下右上）和目标分辨率 `res_deg`，生成严格对齐的经纬度二维网格（`lon_grid`, `lat_grid`），满足 `rasterio` 的 `from_origin` 要求。 | | **`kriging_to_grid()`** | 调用 `pykrige.OrdinaryKriging`，传入点数据和目标网格坐标，执行高斯模型插值。`execute("grid", ...)` 返回与网格形状一致的 2D 数组。 | | **`write_geotiff()`** | 使用 `rasterio` 将插值结果写为标准 GeoTIFF，包含 CRS、transform、压缩（LZW），完全可被 ArcGIS/QGIS 直接加载。 | | **鲁棒性设计** | 自动 CRS 检测与重投影、nodata 处理、异常捕获、进度条、文件名安全化（`p` 替代 `.`）等。 | --- ### ⚠️ 注意事项 & 最佳实践： - **内存优化**：`N_MAX_POINTS=10000` 是平衡精度与速度的经验值。若你的降水场非常平滑（如大范围均匀降水），可降至 `5000`；若地形复杂需高保真，可升至 `20000`，但注意 `pykrige` 对 >20k 点的拟合可能变慢。 - **CRS 必须为地理坐标系**：本代码假设输入为 WGS84（EPSG:4326）。若你的 TIFF 是投影坐标系（如 UTM），**必须先在 ArcGIS/QGIS 中重投影为 WGS84 再运行此脚本**，否则 `0.1°` 无意义。 - **变差函数拟合**：`pykrige` 会自动拟合高斯模型的 `sill`、`range`、`nugget` 参数。如需手动指定（如用前期经验 range），可改用 `variogram_parameters={'sill': ..., 'range': ..., 'nugget': ...}`。 - **并行加速（进阶）**：如需处理数百个文件，可将主循环改为 `concurrent.futures.ProcessPoolExecutor` 并行，但注意 `pykrige` 在多进程下需 `if __name__ == '__main__':` 保护。 ---