Python如何用最近邻法实现精细网格到粗分辨率网格的重采样聚合

DevOps 小满

2026-8-29

问题根因

scipy.interpolate.griddata返回全NaN是必然结果：活跃火点是极稀疏的离散观测，原始1km网格中99%以上的格点为NaN，将包含大量NaN的全量网格传入基于三角网的插值函数时，插值范围无法覆盖稀疏有效点之外的区域，自然全部返回NaN。这类“稀疏点聚合到粗网格”的场景本身就不适合用连续插值方法。
你已写的代码存在两个致命问题：
1. 分开匹配最近纬度、最近经度再拼接索引的逻辑完全错误：经纬度是成对的二维坐标，分开匹配得到的点并非球面上距离火点最近的WRF格点，加上你的MODIS经度网格存在卫星扫描导致的弯曲变形，这种匹配方式会产生公里级的位置偏移。
2. 用Python原生循环逐点查找最近邻效率极低，面对72个时次的近2300万个格点，运行耗时会达到数小时级别，且你没有实现匹配后同格点值聚合求平均的核心步骤。

修正实现方案

采用KD树空间索引批量完成最近邻匹配，再按格点分组聚合求平均，完全匹配需求：仅将有效火点分配至最近的WRF粗格点，无火点区域自动保留NaN，不会丢失任何有效信息，运行效率比原生循环高两个数量级以上。

注：如果WRF输出的lat_wrf/lon_wrf是一维经纬数组，提前用np.meshgrid生成二维网格即可直接运行代码。

import numpy as np
from scipy.spatial import cKDTree

def resample_fire_to_wrf(fine_data, fine_lat2d, fine_lon2d, coarse_lat2d, coarse_lon2d):
    """
    将1km MODIS火点数据聚合降采样到WRF粗网格
    参数:
        fine_data: 单时次精细网格火点数组，维度(4797,4797)，无火点为NaN
        fine_lat2d: 精细网格二维纬度数组，维度(4797,4797)
        fine_lon2d: 精细网格二维经度数组，维度(4797,4797)
        coarse_lat2d: WRF粗网格二维纬度数组，维度(129,109)
        coarse_lon2d: WRF粗网格二维经度数组，维度(129,109)
    返回:
        resampled_data: 重采样后的WRF网格火点数组，维度(129,109)，无火点为NaN
    """
    # 1. 预处理WRF粗网格，构建KD树空间索引
    coarse_rows, coarse_cols = coarse_lat2d.shape
    # 经纬度转弧度，转换为球面三维直角坐标，距离计算等价于大圆距离，避免高纬度经度误差
    coarse_lat_rad = np.deg2rad(coarse_lat2d.ravel())
    coarse_lon_rad = np.deg2rad(coarse_lon2d.ravel())
    coarse_coords = np.vstack([
        np.cos(coarse_lat_rad)*np.cos(coarse_lon_rad),
        np.cos(coarse_lat_rad)*np.sin(coarse_lon_rad),
        np.sin(coarse_lat_rad)
    ]).T
    coarse_tree = cKDTree(coarse_coords)

    # 2. 提取精细网格中所有有效火点
    fire_mask = ~np.isnan(fine_data)
    if not fire_mask.any(): # 当前时次无火点直接返回全NaN数组
        return np.full((coarse_rows, coarse_cols), np.nan)
    fire_vals = fine_data[fire_mask]
    fire_lat_rad = np.deg2rad(fine_lat2d[fire_mask])
    fire_lon_rad = np.deg2rad(fine_lon2d[fire_mask])
    fire_coords = np.vstack([
        np.cos(fire_lat_rad)*np.cos(fire_lon_rad),
        np.cos(fire_lat_rad)*np.sin(fire_lon_rad),
        np.sin(fire_lat_rad)
    ]).T

    # 3. 批量查询每个火点最近的WRF格点索引
    _, nearest_coarse_idx = coarse_tree.query(fire_coords, k=1)

    # 4. 按WRF格点分组求平均
    sum_vals = np.bincount(nearest_coarse_idx, weights=fire_vals, minlength=coarse_rows*coarse_cols)
    count_vals = np.bincount(nearest_coarse_idx, minlength=coarse_rows*coarse_cols)
    # 无火点位置设为NaN
    resampled_flat = np.where(count_vals > 0, sum_vals / count_vals, np.nan)
    # 重构为WRF二维维度
    return resampled_flat.reshape(coarse_rows, coarse_cols)

# ---------------------- 批量处理72个时次数据 ----------------------
# 初始化结果数组，维度(72,129,109)对应[时间×纬度×经度]
wrf_fire_result = np.full((72, 129, 109), np.nan)
for t in range(72):
    wrf_fire_result[t] = resample_fire_to_wrf(
        fine_data=all_time_fire_data[t], # 替换为你的72时次MODIS火点总数组名
        fine_lat2d=lat_mosaic,
        fine_lon2d=lon_mosaic,
        coarse_lat2d=lat_wrf,
        coarse_lon2d=lon_wrf
    )

方案说明

距离计算采用球面三维直角坐标匹配，直接定位球面上最近的格点，避免了单独匹配经纬度的误差，也修正了高纬度经度间隔对应的实际距离缩短的问题，匹配精度满足科研需求。
全程采用numpy向量化操作和cKDTree空间索引，单时次处理通常在1秒内完成，72个时次总耗时不超过2分钟，远快于原生循环。
聚合逻辑完全符合需求：同一WRF格点覆盖范围内的所有火点取平均值，无火点的WRF格点保留NaN，不会出现有效值丢失的问题。

内容的提问来源于stack exchange，提问作者Abhinav Sharma