GrADS第6章变量和函数

1、第第6 6章章 变量和函数变量和函数 在在GrADS运行环境中，参与操作运行环境中，参与操作的数据对象可以是变量、函数和的数据对象可以是变量、函数和表达式形式。表达式形式。 6.1 6.1 变量名变量名 n完全的变量名形式完全的变量名形式n abbrev.file#（dimexpr， dimexpr，）n 其中，其中，abbrev：是数据描述文件中给：是数据描述文件中给出的变量名缩写；出的变量名缩写；file#：为包含此变量的：为包含此变量的已打开的文件序号。已打开的文件序号。 例如：ga-open d:model.ctl open d:model.le.ctl 上例中，上例中，model.c

2、tl文件的序列号为文件的序列号为1，model.le.ctl文件的序列号为文件的序列号为2。序列号。序列号1为缺为缺省情况，可以省略。省情况，可以省略。 括号内的括号内的dimexpr：是对当前维数环境：是对当前维数环境进行的局域维数设置表达式，该表达式仅进行的局域维数设置表达式，该表达式仅对该变量的维数环境进行局域修正，因而对该变量的维数环境进行局域修正，因而仅对固定的维有效。仅对固定的维有效。 设置形式可以采用以下两种方式：设置形式可以采用以下两种方式： 绝对维数表达式为：绝对维数表达式为： X|Y|Z|T|LON|LAT|LEV|TIME = value （value为绝对维数值）；为绝

3、对维数值）；相对维数表达式为：相对维数表达式为： X|Y|Z|T|LON|LAT|LEV|TIME +/- offset （offset为相对于当前维数环境设置的偏差维为相对于当前维数环境设置的偏差维数值）。数值）。 例如例如： z.3（lev=500）tv.1（time-12hr）rh 缺省的当前文件中的变量缺省的当前文件中的变量rhq.2（t-1,lev=850）2号文件中相对于当前时刻的号文件中相对于当前时刻的前一时刻，高度为前一时刻，高度为850hPa等等压面上的变量压面上的变量q表示文件表示文件3中高度为中高度为 500hPa 等压面上的变量等压面上的变量Z相对于当前时刻之前相对于当

4、前时刻之前12小时小时时刻的时刻的1号文件中的变量号文件中的变量tvndset model.le.dat noptions little_endian cray_32bit_ieeenUNDEF -2.56E33 nTITLE 5 Days of Sample Model OutputnXDEF 72 LINEAR 0.0 5.0 nYDEF 46 LINEAR -90.0 4.0 nZDEF 7 LEVELS 1000 850 700 500 300 200 100nTDEF 5 LINEAR 02JAN1987 1DY nvars 8 nps 0 99 Surface Pressurenu

5、 7 99 U Windsnv 7 99 V Windsnz 7 99 Geopotential Heightsnt 7 99 Temperaturenq 5 99 Specific Humiditynts 0 99 Surface Temperaturenp 0 99 PrecipitationnENDVARS 2. 举例说明举例说明ndset model.grbntitle Sample Model Data for lats4d Tutorialnundef 1e+20ndtype gribnindex model.gmpnxdef 72 linear 0.000000 5.000000

6、nydef 46 linear -90.000000 4.000000nzdef 7 levels 1000 850 700 500 300 200 100 ntdef 5 linear 0Z1jan1987 1dynvars 8nps 0 1, 1, 0, 0 Surface pressure hPanua 7 33,100 Eastward wind m/snva 7 34,100 Northward wind m/snzg 7 7,100 Geopotential height mnta 7 11,100 Air Temperature Knhus 7 51,100 Specific h

7、umidity kg/kgnts 0 11,105, 2 Surface (2m) air temperature Knpr 0 59, 1, 0, 0 Total precipitation rate kg/(m2*s)nendvarsga-open d:model.ctl open d:model.le.ctl set t 1 将时间维数固定在第将时间维数固定在第1个时次个时次 d z.2(lev=500) set z 3 设置垂直高度层次为第设置垂直高度层次为第3个层次个层次,即即700hPa d z(z+1)同时打开两个文件同时打开两个文件: model.ctl文件序列号为文件序列号为

8、1，model.le.ctl 文件序列号为文件序列号为2.显示显示2号文件中号文件中500hPa等压面等压面 上的变量上的变量z(绝对维数绝对维数)显示显示1号文件中相对于第号文件中相对于第3个层次的后一个层次的后一个层次，即个层次，即500hPa面上的变量面上的变量z （相对（相对维数）维数）6.2 6.2 表达式表达式1 1、表达式的构成、表达式的构成n运算符：运算符： + + 加，加，- - 减，减，* * 乘，乘，/ / 除；除；n运算域：运算域： 变量变量, , 函数和常数；函数和常数；n括号：用于控制运算的次序。括号：用于控制运算的次序。2. 2. 运算规则运算规则n对相同网格点上

9、的不同变量进行运算时，对相同网格点上的不同变量进行运算时，只要有一个变量在某格点的值为缺测，则只要有一个变量在某格点的值为缺测，则该网格点的运算结果为缺测值；当被该网格点的运算结果为缺测值；当被0 0除除时的结果也设为缺测。时的结果也设为缺测。n对多个数据文件的数据做运算时，两种数对多个数据文件的数据做运算时，两种数据网格要一致，即运算对象的维数变化范据网格要一致，即运算对象的维数变化范围要一致。围要一致。n如某一格点的变动维数多于其他格点，则如某一格点的变动维数多于其他格点，则具有较少变动维数的格点其维数环境将被具有较少变动维数的格点其维数环境将被拓展，以便于运算操作。拓展，以便于运算操作。

10、例如：例如：z - z(t-1)t(lev=200)-t(lev=850) ave(z,t=1,t=5) z - ave(z,lon=0,lon=360,-b)变量变量Z的时间变化的时间变化200hPa与与850hPa等压面上等压面上的变量的变量t的变化的变化变量变量Z从第从第1时次到第时次到第5时次的平均时次的平均变量变量Z 的纬向偏差的纬向偏差6.3 6.3 函数函数1、调用方式、调用方式 采用直接引用方式，形式如：采用直接引用方式，形式如： 函数名函数名(参数参数1,参数参数2,)2. 2. 常用函数常用函数（1）ave函数函数n格式：格式： ave(expr,dexpr1,dexpr2

11、,tincr)n功能：通用的求平均函数。功能：通用的求平均函数。n说明：说明：expr是由是由dexpr1和和dexpr2定义的维数定义的维数 范围内求平均的量。范围内求平均的量。 例如：ga-open model.ctl set lev 500 d ave(z,t=1,t=5) 显示显示500hPa等压面上从第等压面上从第1时次到第时次到第5时次时次变量变量z的平均量的平均量set parea 1 8 1 6set grid offset grads offset xlint 20set xlopts 1 2 0.18set ylopts 1 2 0.18set lon 60 180set

12、lat 0 60set lev 200set t 1define a=ave(v,t=7,t=660,12)d aprintc例如：例如： 求纬向平均：求纬向平均： 如果如果ave(expr,lon=0,lon=360)在端点求在端点求 了两次平均；了两次平均； 如若使终端点取半数权重，可用：如若使终端点取半数权重，可用： ave(expr,lon=0,lon=360,-b)。 set xyrev onset lev 500set t 1set lat 0 90 set x 1define a=ave(v,lon=0,lon=360,-b)d aset cmark 0set ccolor 2d

13、 const(a,0)(2) const函数函数n格式：格式：const（expr,constant）n功能：设置部分网格点的值取为常数功能：设置部分网格点的值取为常数constantn说明：所有非缺测格点处的说明：所有非缺测格点处的expr值取为常数，值取为常数，flag为选项为选项，如果加上选项，如果加上选项-a，则所有网格点，则所有网格点值均设定为指定的常数，值均设定为指定的常数，如果加上选项如果加上选项-u，则，则只把只把缺测格点处缺测格点处的的expr值设定为常数。值设定为常数。n该函数对格点和台站资料均适用。该函数对格点和台站资料均适用。（3）aave函数函数n格式：格式： aav

14、e(expr,xdim1,xdim2,ydim1,ydim2)n功能：求面积平均。功能：求面积平均。n说明：说明：expr: 任何表达式任何表达式xdim1: X 维数方向的起始维数表达式维数方向的起始维数表达式xdim2: X 维数方向的终止维数表达式维数方向的终止维数表达式ydim1: Y 维数方向的起始维数表达式维数方向的起始维数表达式ydim2: Y 维数方向的终止维数表达式维数方向的终止维数表达式set lev 850set t 1 48set x 1set y 1define a=aave(v,lon=100,lon=120,lat=30,lat=40)d a（4）mag函数函数n

15、格式：格式：mag(aexpr,bexpr)n功能：计算功能：计算n说明：说明：aexpr表示表示a变量，变量，bexpr表示表示b变量，变量，该函数对格点和台站资料均适合。该函数对格点和台站资料均适合。22)exp()exp(rbraset lon 60 200set lat -20 80d mag(u.1,v.2)（5）hcurl函数函数n格式：格式：hcurl(uexpr,vexpr)n功能：计算功能：计算垂直涡度垂直涡度。n说明：说明： uexpr表示表示U风速分量，风速分量，vexpr表示表示V风速分风速分量，风速单位用量，风速单位用m/s。边界上的涡度值设定。边界上的涡度值设定为缺

16、测。为缺测。set lon 60 200set lat -20 80 d hcurl(u,v) （6）hdivg函数函数n格式：格式：hdivg(uexpr,vexpr)n功能：计算功能：计算水平散度水平散度。n说明：说明： uexpr表示表示U风速分量，风速分量，vexpr表示表示V风速分风速分量，风速单位用量，风速单位用m/s。set lon 60 200set lat -20 80 d hdivg (u,v) （7）skip函数函数n格式：格式：skip（expr,skipx,skipy）n功能：设定样本的取样密度。功能：设定样本的取样密度。n说明：说明：skipx,skipy 数值决定

17、数值决定X和和Y方向的取方向的取样密度（取值样密度（取值1可以省略不给）可以省略不给）*该函数主要用于对该函数主要用于对矢量场的稀疏化矢量场的稀疏化显示。显示。 d u.1;v.2d skip(u.1,2,2);skip(v.2,2,2)（8）cdiff函数n格式：格式：cdiff（expr,dim）n功能：中央差分函数n说明：expr表示进行差分运算的量， dim表示进行差分运算的维数方向， 为X，Y，Z中的任一个字符。 边界格点的运算值设为缺测。边界格点的运算值设为缺测。n例如：计算温度平流define dtx=cdiff(t,x)define dty=cdiff(t,y)define d

18、x=cdiff(lon,x)*3.1416/180define dy=cdiff(lat,y)*3.1416/180d -1*(u*dtx)/(cos(lat*3.1416/180)*dx)+v*dty/dy)/6.37e6（9 9）smth9smth9函数函数格式：格式：smth9(expr)smth9(expr)功能：功能：9 9点平滑函数点平滑函数说明：说明：exprexpr表示须平滑的表达式。表示须平滑的表达式。 该函数主要用于改善图形的输出质量，该函数主要用于改善图形的输出质量， 使其更加美观。使其更加美观。当格点资料为当格点资料为1 1维时，维时，9 9点平滑退化为点平滑退化为3

19、3点平滑。点平滑。d smth9(b1)（11）vint函数n格式：vint(psexpr,expr,top)n功能：质量加权垂直积分函数n说明：psexpr表示地面气压变量， expr表示被积变量， top表示积分上限(须为常数，hPa) 本函数要求定义水平维数环境。Vint( psexpr, expr, top)=f/g* sum(expr * Delta(level)f=100，g=9.8，sum是求和的意思是求和的意思 n例如：计算整层水汽通量fuq;fvq。open e:bsuq.ctlopen e:bspres.ctlset gxout fwriteset fwrite e:bsm

20、onzuq.grdi=1while(i=660)set t i set x 1 144set y 1 73set z 1define fuq=vint(v.2(z=1)/100,v,300)set t iset x 1 144set y 1 73set z 1d fuqi=i+1 endwhiledisable fwritereinit1998年夏季整层水汽通量矢量（箭头）、水汽通量值（阴影）和水汽路径（黑色曲线）站点画图站点画图（1）写成）写成站点格式站点格式grd数据；数据；（2）生成）生成map文件文件；（3）建立）建立格点文件格点文件（建立后可以继续使用）；（建立后可以继续使用）；（4

21、）将站点资料通过）将站点资料通过插值函数插值函数插值到格点文件插值到格点文件上，然后再画图。上，然后再画图。 文件的时间一定要与站点数据时间一致文件的时间一定要与站点数据时间一致。 范围：范围：lat：15 55 N lon：70 140 E 格距：格距：1*1ngrid.grd文件的每个点上均赋值为文件的每个点上均赋值为1格点文件grid.grd ：第二章数据处理中第二章数据处理中70-73页写出的格点数据页写出的格点数据grid.grd，对应对应ctl为为grid.ctl与该格点文件相对应的数据描述文件grid.ctl如下：dset d:datagrid.grdundef -999.0ti

22、tle Sample GRID Dataxdef 71 linear 70 1ydef 41 linear 15 1zdef 1 linear 500 1tdef 1 linear jul1951 1movars 1g 0 99 grid data endvars有了这个格点文件及与其相对应的数据描述文件后，就可以通过有了这个格点文件及与其相对应的数据描述文件后，就可以通过GrADS插值函数将站点资料插值到该格点文件上然后使用了。插值函数将站点资料插值到该格点文件上然后使用了。 (12) oacres函数函数n格式：格式：oacres(gexpr, sexpr ) n功能：功能： cressm

23、an 插值函数插值函数n说明：说明： 常用于常用于 站点资料插值为格点资料。站点资料插值为格点资料。 gexpr 大于站点数据范围的格点数据变量名大于站点数据范围的格点数据变量名sexpr 站点数据变量名站点数据变量名radii 影响半径，影响半径， 默认值为默认值为: 10,7,4,2,1 n为避免插值出现虚假结果为避免插值出现虚假结果, , 参考网格的间参考网格的间隔与站点间距相近为好隔与站点间距相近为好将站点资料通过插值函数将站点资料通过插值函数插值插值到到某个格点文件某个格点文件上，上，然后再使用。然后再使用。与该格点文件相对应的数据描述文件grid.ctl如下：dset d:data

24、grid.grdundef -999.0title Sample GRID Dataxdef 71 linear 70 1ydef 41 linear 15 1zdef 1 linear 500 1tdef 1 linear jul1951 1movars 1g 0 99 grid data endvarsdefine a=oacres(g,rain.2) 站点数据变量名站点数据变量名open d:datagrid.ctl open d:datarain7.ctl （13）maskout函数函数n格式：格式：maskout(expr,mask) n功能：功能： 标记函数。当标记函数。当 mas

25、k所在的网格点值取为所在的网格点值取为负值时，对应格点上的负值时，对应格点上的 expr的值设为缺测值，不的值设为缺测值，不参与运算或画图。参与运算或画图。 即：只画出即：只画出mark大于大于0的的expr。n说明：说明： expr 需要处理的变量名需要处理的变量名mask 标记变量名标记变量名 常用于输出常用于输出 mask 代表的陆地区或海洋区；代表的陆地区或海洋区；也常用于我国站点资料绘图。也常用于我国站点资料绘图。open d:datagrid.ctlopen d:datar7.qht.ctlset grads offenable print d:datarain7.qht.gmfs

26、et mpdset hires cnworldset lon 73 135set lat 15 55define a=oacres(g(t=1),p.2,1.5) （“oacres”为插值函数）为插值函数）define a1=maskout(a,g-0.5) （“maskout”为标记函数）为标记函数）define aa=smth9(a1) （“smth9”9点平滑函数）点平滑函数）set grid offset gxout shadedset black -0.2976 0.2976d aaset gxout contourset cint 0.1d aa一般将格点数据作为第一个打开的文一般

27、将格点数据作为第一个打开的文件，在当前件，在当前缺省文件为格点数据缺省文件为格点数据时时, 可可以用以用define命令将站点数据的客观分命令将站点数据的客观分析结果保存到变量中。析结果保存到变量中。printdisable printreinit;open e:datagrid.ctlopen e:datat.ctlset lon 70 140set lat 15 55define aa=oacres(g.1,t1.2,50,40,30,20,10,5,2)define aaa=maskout(aa,g.1-0.5)set gxout shadedd aaaset gxout contour

28、d aaaprintim e:datasample.png whitereinit；n地形高度资料：dxgd.ctl 变量名：hn700h Pa上，地形大于3000米以上的区域是没有风场的，所以可以利用maskout函数将该地区风场不绘制出来reinitopen E:bsmonu.ctlopen E:bsmonv.ctlopen e:dxgd.ctlenable print E:uv.gmfset lev 700set t 1set lon 70 140set lat 10 60d maskout(v,3000-h.3(z=1,t=1);maskout(v.2,3000-h.3(z=1,t=1

29、)printdisable printreinit;dset e:dxgd.grdundef -9.99E+33xdef 144 linear 0 2.5ydef 73 linear -90 2.5zdef 1 linear 1 1tdef 1 linear 01jan1951 1movars 1h 1 99 hendvarsset t 1set lon 70 140set lat 10 60set gxout shadedset clevs 3000set ccols 0 1d h.3(z=1)set lev 700set t 1set lon 70 140set lat 10 60set

30、gxout vectord v;v.2printdisable printreinit;其他函数n数学函数数学函数 统计函数：统计函数： 网格处理函数n常数 const (expr, const, | -u | -a )n剔除 maskout (expr, mask_expr )n跳点skip (expr, skip_x, skip_y )n加权9点平滑点平滑 smth9 ( expr ) 1-D 3点平滑点平滑n中央差分 cdiff ( expr, dim )物理量计算函数虚温,相对湿度求比湿相对湿度求比湿-虚温,相对湿度求温度相对湿度求温度-垂直涡度-中央差求水平散度-全风速-6.4 6.

31、4 定义变量定义变量1、定义形式、定义形式define varname = expr 其中，其中，varname为新变量名，为新变量名，expr为表达为表达式。所定义的新变量可以用于随后的表达式。所定义的新变量可以用于随后的表达式中。式中。define命令命令只只适用于适用于格点数据格点数据2. 存储形式存储形式 新定义的变量新定义的变量varname只存在于内存中，只存在于内存中，所以建议不要定义过大的维数范围。所以建议不要定义过大的维数范围。 3.3.维数环境的设置维数环境的设置n用户定义的变量可以有用户定义的变量可以有04个变化的维数。个变化的维数。假设定义一个假设定义一个4维的变量，在

32、维的变量，在“define”命命令执行后，记住改变维数环境，使之变化令执行后，记住改变维数环境，使之变化的维数小于的维数小于4。n当当4个个维数同时变化时，维数同时变化时，“define”命令是命令是GrADS中所有命令中中所有命令中唯一唯一可行的命令。命令。 open e:uv.ctlset lon 60 150set lat 0 40set lev 850 200set t 1 48define temp=mag(u,v)set t 5set lev 200d tempcreinit定义后使用时改变其维数环境定义后使用时改变其维数环境定义数组维数环境（定义数组维数环境（4维）维）此时显示的

33、变量此时显示的变量temp，其时间和层次维是固定的，其时间和层次维是固定的定义数组定义数组n如果定义的变量具有某些固定的维数，随如果定义的变量具有某些固定的维数，随后使用这个变量，无论维数环境如何改变，后使用这个变量，无论维数环境如何改变，此变量在固定维数上保持不变。此变量在固定维数上保持不变。 set lon -180 0 set lat 0 90 set lev 500 set t 1 define zave = ave(z,t=1,t=5) X、Y维度不固定Z、t维度固定set t 5 set lev 200 d zave 显示的是显示的是500 h Pa，1-5时刻平均时刻平均的值的值

34、n当定义的变量具有某些变化的维数环当定义的变量具有某些变化的维数环境，然后再将这些维数固定起来，这境，然后再将这些维数固定起来，这时所得的变量将被固定在该维数上。时所得的变量将被固定在该维数上。set lon 0 180set lat -90 90 set lev 500 set t 5 define temp = z set lat 40 d temp Y维固定为40N，显示的只能是40N 的值set lon 200set lat 40 d temp 超出定义范围，超出定义范围，显示为缺测显示为缺测n用户也可以使用局地维数环境。用户也可以使用局地维数环境。n如果局地环境是在定义变量时变化的如果局地环境是在定义变量时变化的维数维数环境内环境内，则，则显示显示局地维数环境下的局地维数环境下的变量值变量值。n如果这个维数在定义变量时是如果这个维数在定义变量时是固定的固定