Python实现纬度加权平均

hhhz

XiaoMaFenJu 发表于 2024-4-22 20:58
开方肯定是有影响的，但是我也不知道为什么会开方（我也见过有开方的）

看到有做EOF之前纬度加权开方，是因为EOF分析时会对X方向和Y方向同时进行分解，所以XY方向都会会进行加权，两个方向加权后的面积权重时cos纬度，单个方向就需要开根号https://bbs.06climate.com/forum.php?mod=viewthread&tid=36756

hhhz

XiaoMaFenJu 发表于 2024-4-22 20:58
开方肯定是有影响的，但是我也不知道为什么会开方（我也见过有开方的）

找到有EOF对数据加权用了开方，这个帖子下解释了原因：https://bbs.06climate.com/forum.php?mod=viewthread&tid=36756      因为EOF分析时会对X方向和Y方向同时进行分解，所以XY方向都会会进行加权，两个方向加权后的面积权重时cos纬度，单个方向就需要开根号

XiaoMaFenJu

hhhz 发表于 2024-4-24 21:40
找到有EOF对数据加权用了开方，这个帖子下解释了原因：https://bbs.06climate.com/forum.php?mod=viewthr ...

那如果不是eof也可以直接cos？不过开不开根的误差应该都不会很大

hhhz

是的，可能开方是针对某些特定数据处理要做的吧

姚小蘭

参考楼主的代码，三维的我这么处理：

def AreaWeightMean_3D(data3D, lat, lon):
    '''
    data3D: 要进行区域加权平均的变量  3D: [time, lat, lon]
    lat: data3D对应的纬度 1D
    lon: data3D对应的经度 1D
    '''
    jlat = lat.shape[0]
    rad = 4.0 * np.arctan(1.0) / 180.0
    re = 6371220.0
    rr = re * rad
   
    # 经度间隔，因为是均匀分布，可以直接使用第一个间隔作为所有间隔
    dlon = np.abs(lon[2] - lon[1]) * rr
    dx = dlon * np.cos(lat * rad)
   
    # 纬度间隔不均匀，需要单独计算每个点的间隔
    dy = np.zeros(jlat)
    dy[0] = np.abs(lat[2] - lat[1]) * rr
    dy[1: jlat - 1]  = np.abs(lat[2: jlat].values-lat[0: jlat - 2].values)*rr * 0.5
    dy[jlat - 1] = abs(lat[jlat - 1] - lat[jlat - 2]) * rr
    area = dx * dy
        
    # 扩展 area 维度以匹配 data3D
    area_3D = np.expand_dims(area, axis=(0, 2))

    # 处理缺失值
    valid_mask = ~np.isnan(data3D)
    data3D = np.nan_to_num(data3D)
    # 计算加权和
    sumtop = np.sum(data3D * area_3D * valid_mask, axis=(1, 2))
    sumbottom = np.sum(area_3D * valid_mask, axis=(1, 2))

    dataAreaWeightMean = sumtop / sumbottom
   
    return dataAreaWeightMean