Python实现纬度加权平均

XiaoMaFenJu

原代码来自：

https://blog.csdn.net/mihuanchengbao/article/details/127572695

因python中没有现成的包去计算面积加权平均（即NCL中的wgt_areaave），故找到作者"可乐要加冰_ice"的代码，但计算时间过长，故在此利用numpy进行改进，代码如下：

1. def AreaWeightMean(data2D, lat, lon):
   
2.     '''
   
3.     data2D: 要进行区域加权平均的变量  2D: [lat, lon]
   
4.     lat: data2D对应的纬度 1D -90° 和 90°  因为NCL 和 Python 计算 np.cos(90 * rad) 值差的很大，故使用89.9°进行替代计算
   
5.     lon: data2D对应的经度 1D
   
6.     '''
   
7.     if lat.data[0] == 90:
   
8.         lat.data[0] = 89.9
   
9.         pass
   
10.     if lat.data[-1] == -90:
    
11.         lat.data[-1] = -89.9
    
12.         pass
    
13. 
    
14.     jlat = lat.shape[0]
    
15.     rad = 4.0 * np.arctan(1.0) / 180.0
    
16.     re = 6371220.0
    
17.     rr = re * rad
    
18.     dlon = np.abs(lon[1] - lon[0]) * rr
    
19.     dx = dlon * np.cos(lat * rad)
    
20.     dy = np.zeros(jlat)
    
21.     dy[0] = np.abs(lat[1] - lat[0]) * rr
    
22.     dy[1: jlat - 1] = np.abs(lat[2: jlat].values-lat[0: jlat - 2].values)*rr * 0.5
    
23.     dy[jlat - 1] = abs(lat[jlat - 1] - lat[jlat - 2]) * rr
    
24.     area = dx * dy
    
25.    
    
26.     weight = np.expand_dims(area, 1).repeat(len(lon), axis=1)
    
27.     new_data = np.average(data2D, weights=weight)
    
28.     return new_data

复制代码

本代码改进地方有二：一是在纬度存在90°/-90°时，因极点处权重差距不合理，选取89.9°进行替代计算；二是利用np.average()函数加快计算过程；三是源代码在dy[0]处有问题，已更改


=======================2023.11.04============================
面积加权本质就是随纬度升高、平面半径有所减小，重点就在cos(lat)上，故进行重写并忽略无关因素，且支持3维数据，代码如下：

1. def Area_Mean(data, lat, lon):
   
2.     '''
   
3.     by XiaoMaFenJu
   
4.     data: 要进行区域加权平均的变量，支持2、3维  2D: [lat, lon]  3D：[time, lat, lon]
   
5.     lat: data2D对应的纬度 1D
   
6.     lon: data2D对应的经度 1D
   
7.     '''
   
8. 
   
9.     if data.ndim == 2:
   
10.         y_weight2D = abs(np.cos(lat*np.pi/180))
    
11.         weight2D = np.expand_dims(y_weight2D, 1).repeat(len(lon), axis=1)
    
12.         # print(weight2D)
    
13.         new_data = np.average(data, weights=weight2D)
    
14.         return new_data
    
15.     elif data.ndim == 3:
    
16.         y_weight2D = abs(np.cos(lat*np.pi/180))
    
17.         weight2D = np.expand_dims(y_weight2D, 1).repeat(len(lon), axis=1)
    
18.         weight3D = np.expand_dims(weight2D,0).repeat(len(data[:,0,0]),axis=0)
    
19.         print(weight3D.shape)
    
20.         new_data = np.average(data, weights=weight3D, axis=(-1, -2))
    
21.         return new_data
    
22.     else:
    
23.         print('输入数据非2&3维')

复制代码

新函数计算所得结果与老函数几乎相同。对于全球气温数据来说90°的差异造成的影响很小，故新函数就没有对90°进行特别处理。

sunyuqing

krwlng 发表于 2023-11-11 09:48
想问一下楼主如果数据有缺测的化需要怎么预处理呢？

可以在循环里加：mask = ~np.isnan(data)      ; data_mean = np.average(data[mask], weights=wgts[mask])  大概这样

XiaoMaFenJu

顶一下

krwlng

想问一下楼主如果数据有缺测的化需要怎么预处理呢？

XiaoMaFenJu

krwlng 发表于 2023-11-11 09:48
想问一下楼主如果数据有缺测的化需要怎么预处理呢？

我暂时没有这个需求呢，我觉得应该就是np.average跟nan的问题，你可以参考https://www.coder.work/article/371500

krwlng

XiaoMaFenJu 发表于 2023-11-11 12:45
我暂时没有这个需求呢，我觉得应该就是np.average跟nan的问题，你可以参考https://www.coder.work/articl ...

感谢回复！

krwlng

sunyuqing 发表于 2023-11-12 22:28
可以在循环里加：mask = ~np.isnan(data)      ; data_mean = np.average(data[mask], weights=wgts[mask ...

感谢帮助！我加了这句画出我想要的图了，非常感谢！

lilu92

学习了！收藏验证一下

hhhz

我看到有的在做纬度加权时候对纬度cos进行了开方作为每个纬度的权重，请问这个开不开方有什么影响嘛

XiaoMaFenJu

hhhz 发表于 2024-4-22 11:09
我看到有的在做纬度加权时候对纬度cos进行了开方作为每个纬度的权重，请问这个开不开方有什么影响嘛

开方肯定是有影响的，但是我也不知道为什么会开方（我也见过有开方的）