赞
踩
当处理全球数据的时候,最常见的分析就是全球、北半球平均了,这时候一定要注意要根据纬度等信息进行加权,求加权平均,要不然结果肯定会不准。
对于NCL中,有现成的函数实现:wgt_areaave
但是对于Python中,常规的计算包里面是没有的,或许现在有其他包可以实现,也可以告诉我~。于是基于numpy写了下面的函数:
- def AreaWeightMean(data2D, lat, lon):
- '''
- data2D: 要进行区域加权平均的变量 2D: [lat, lon]
- lat: data2D对应的纬度 1D 最好不要包括-90° 和 90° 因为NCL 和 Python 计算 np.cos(90 * rad) 值差的很大
- lon: data2D对应的经度 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]-lat[0: jlat - 2])*rr * 0.5
- dy[jlat - 1] = abs(lat[jlat - 1] - lat[jlat - 2]) * rr
- area = dx * dy
- # dataAreaWeightMean = np.sum(np.dot(area, data2D))/np.sum(area)
- sumtop = 0
- sumbottom = 0
- for id1 in range(data2D.shape[0]):
- for id2 in range(data2D.shape[1]):
- sumtop = sumtop + data2D[id1, id2]*area[id1]
- sumbottom = sumbottom +area[id1]
- dataAreaWeightMean = sumtop/sumbottom
- return dataAreaWeightMean
'运行
看一下结果对比吧:
直接不加权平均(np.mean(data3D, axis = (1,2)))得到的结果:
加权平均之后的:
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。