关于grads加权平均/求和的问题

http://disc.sci.gsfc.nasa.gov/giovanni/giovanni_analysis_functions.shtml

加权平均ave(expr, dim_expr1,dim_expr2, tincr, { |-b } )

加权区域平均aave(expr, { xdim1,xdim2,ydim1,ydim2 | global | g } )

加权求和（一维求和）sum(expr, dim_expr1,dim_expr2, tincr, { |-b } )（同ave）

加权区域和（二位求和） asum (expr, { xdim1,xdim2, ydim1, ydim2| global | g } )

ave(expr, dim1, dim2<,tinc> <,-b>)

Usage Notes

1平均所用的网格坐标来自于缺省文件。若dim1和dim2为地理坐标（即经纬度），则按缺省文件定义的网格，将地理坐标转换到距离最近的整数网格坐标。

2若没有给出边界标识–b，则始末点的权重与内部点一样，即始末点实际是离dim1和dim2最近的网格坐标点，若给出边界标识-b，则始末点为给定值。

3对于非线性的网格间距，此平均函数以网格间距为权重系数。对应沿经线方向上的平均（不同纬度），权重系数为网格南北边界纬度的正弦函数差，而网格边界通常定义为两相邻网格点的中点。

Example

1若求纬向平均（沿纬线方向）ave(z,lon=0,lon=360)

则地理坐标将被转换成网格坐标，在这里X实际取1到181，0°（和360°）在平均中被使用了两次，为了使平均过程中的始末点有更合适的权重，可使用-b，ave(z,lon=0,lon=360,-b)或者直接使用网格坐标ave(z,x=1,x=180)

ave(z,lon=0,lon=360)平均的前后边界为-1→361（359-1）

ave(z,lon=0,lon=360,-b)平均的前后边界为0→360

2平均算子可如下嵌套ave(ave(z,x=1,x=180),y=1,y=46)这个嵌套可以用来求区域平均。但对于区域平均，aave更实用。

当将平均算子进行嵌套时，嵌套的顺序对结果有很大的影响。在grads中，X维变化最快，其次是Y维，其次Z维，最后是T维。当进行平均嵌套时，将快变维放在最里层。如ave(ave(t,x=1,x=180),t=1,t=20)比ave(ave(t,t=1,t=20),x=1,x=180)更高效。

3使用define命令可以提高特定的算子的效率。如若要计算标准差，算法若写作：

sqrt(ave(pow(ave(z,t=1,t=20)-z,2),t=1,t=20))

虽然也没错，但是里层的平均ave(z,t=1,t=20)-z,2)将会重复计算20次。若将里层的平均预先定义，计算速度将会大幅提高。

define zave =ave(z,t=1,t=20)

dsqrt(ave(pow(zave-z,2),t=1,t=20))

aave(expr, xdim1, xdim2, ydim1, ydim2)

Usage Notes

1对于全球平均，可简化为aave(expr,g)或aave(expr,global)，结果与

aave(expr, lon=0, lon=360, lat=-90,lat=90) 相同。

2 若无缺测值，则aave与X和Y方向两次ave嵌套的结果相同，但aave更快。即：ave(ave(expr,x=1,x=72),y=1,y=46)与aave(expr,x=1,x=72,y=1,y=46)结果相同。

但当有缺测值时，二者结果不同。考虑如下网格，其中U代表缺测。

6      18       3       5

10      10      10     10                                                                      

12       U      U       U

若使用ave嵌套，X维的ave为8,10,12，再进行Y维ave得10

若使用aave，得84/9=9.33

3函数aave默认X方向为经度，Y方向为纬度，在纬度方向有权重，权重因子为格点南北边界的纬度的正弦之差。

4 与ave不同，aave在经纬度坐标下默认使用精确边界进行平均计算。若给定区域的边界值与网格格点边界不重合，则边界值根据格点中心点离给定边界的距离进行计算。

5 若用格点坐标定义区域，则格点坐标将会转换为地理坐标以确定精确边界。这种转换会根据缺省文件的网格定义进行。要注意的是，这种转换使用的是网格边界点而不是中心点。例如：aave(expr,x=1,x=72,y=1,y=46)

在区域定义中x=1,x=72,y=1,y=46分别代表1E，72E，1N，46N，但实际代表区域0.5E-72.5E，0.5N-46.5N。这4个网格边界点将会根据缺省文件中定义的网格转换成经纬度坐标。假设x=1代表0E，x=72代表355E，则区域的边界为-2.5和357.5，这恰好覆盖了全部的经度。在Y方向，当边界超过极点时，aave函数会识别并给以适当权重。

6在asum中没有纬度方向的权重。atot有。

clc;clear;

yc=[0:5:90];

ys=yc-2.5;

yn=yc+2.5;

yc=yc./180*pi;

ys=ys./180*pi;

yn=yn./180*pi;

a=cos(yc);

b=sin(yn)-sin(ys);

s1=a/sum(a);

s2=b/sum(b);

%sin(yn)-sin(ys)=(yn-ys)*cos(y)

syms ya;

limit((sin(y+a/2)-sin(y-a/2))/(a*cos(y)),a,0,'right')