Grads 笔记

注：2012-2016年抄的，已经不用这软件，停止更新

Grads绘图程序几例

例1 （相关及信度检验；中国160站插值绘图）

中国16站插值绘图；打开多个文件；设置画布大小；等值线阴影叠加绘图；改等值线间隔，加粗；

部分数值绘制阴影，自定义颜色；调用长江黄河、中国底图；不打印grads标注；填加标题

'reinit'

'open D:\zother\grid.ctl'

'open D:\zother\R_autumn.ctl'

'open D:\zother\Sig_autumn.ctl'

'enable print D:\zother\R_Sig_autumn.gmf'

*设置画布大小

'set parea 1 10 1 8.4'

*不打印grads标注

'set grads off'

 

*要先画阴影后画等值线，否则会引起覆盖！

*---调用地图---------------

'set mpdset cn cnriver'

'set mpdset cnworld'

'set mpdset cnscs'

*----设置经纬度------------

'set lon 70 135.5'

'set lat 17 55'

*--------------------------

*引用第3号文件中的变量，要在后面加".3"

'define a2=oacres(g,Sig_autumn.3)'

'define b2=maskout(a2,g-0.5)'

'define bb2=smth9(b2)'

'set xlopts 1 10 0.18'

'set ylopts 1 10 0.18'

*-------阴影图--------

'set gxout shaded'

*设置自定义的彩虹色序列：蓝、中蓝、白

'set rbcols 4 11 0'

*只画0.1和0.05的值

'set clevs 0.05 0.1'

'd bb2'

 

*-------画R线型图--------

*设置地图颜色、粗细、线型

'set map 1 1 6'

*引用第2号文件中的变量，要在后面加".2"

'define a1=oacres(g,R_autumn.2)'

*--------------------------

'define b1=maskout(a1,g-0.5)'

'define bb1=smth9(b1)'

'set xlopts 1 10 0.18'

'set ylopts 1 10 0.18'

*--------------------------

*-------线型图--------

*设置颜色、粗细

'set ccolor 1'

'set cthick 8'

'set gxout contour'

*等值线间隔0.1

'set cint 0.1'

'd bb1'

*--------------------------

 

'draw title Corr. (Nino3.4,rainfall) Autumn-Autumn'

*--------------------------

*'cbarn.gs'

*色棒大小 类型 左右 上下

'cbarn 1 0 5.6 0.9'

 

'set cthick 7'

'set clopts 1 6 0.1'

**'d:\zother\southsea.gs'

'print'

'disable print'

;

例2（降水合成图；中国160站插值绘图）

类似例1，为把不同图色棒调一致，画指定等值线

 

'reinit'

'open D:\zother\grid.ctl'

'open D:\zother\Spring_Composite_3.4_above1.ctl'

'enable print D:\zother\Next_Spring_Composite_3.4_above1.gmf'

*不打印grads标注

'set grads off'

'set lon 75 135.5'

'set lat 17 55'

'set mpdset cn cnriver'

'set mpdset cnworld'

'set mpdset cnscs'

'define a1=oacres(g,Sp_Comp.2)'

'define b1=maskout(a1,g-0.5)'

'define bb1=smth9(b1)'

'set xlopts 1 10 0.18'

'set ylopts 1 10 0.18'

'set gxout shaded'

'set clevs -40 -30 -20 -10 0 10 20 30 40 50 60'

'd bb1'

'set gxout contour'

'set clevs -40 -30 -20 -10 0 10 20 30 40 50 60'

'd bb1'

*'draw title  '

'cbarn.gs'

'set cthick 8'

'set clopts 1 6 0.1'

**'d:\zother\southsea.gs'

'print'

'disable print'

;

 

Grads打开及读取nc文件

1.版本要求：太低版本的打不开，可下载气象家园整合版

2.打开命令：sdfopen XXX.nc    ；查看文件ctl详细信息：q ctlinfo

例如：

上图说明：

xdef ：144个格点；起始格点：0；格距：2.5

ydef：73个格点；起始格点：-90；格距：2.5

zdef：17个层次，分贝为：1000,925......

tdef：820个时次；起始时间为：00Z01JAN1948；时间间隔：一个月

vars：变量数：1个

变量名：hgt

存储顺序：t,z,y,x （读写文件时要按此顺序写循环）

 

3.提取所需资料的gs文件

（以上图为例）

'sdfopen E:\ncep\hgt.mon.mean.nc'  （打开的nc文件）

'set gxout fwrite'

'set fwrite E:\ncep\hgt200.grd'    （提取的grd文件名称及储蓄路径）

'set x 1 144'    （根据所打开的文件修改）

'set y 1 73'      （根据所打开的文件修改）

'set lev 200'    （根据所需要提取的层次修改，此处提取200hpa的）

'set t 1 820'     （根据所打开的文件/所需修改）

'd hgt'               （根据所打开的文件/所需修改，变量名见ctl）

'disable fwrite'

;

 

 

Grads配ctl文件画图

常用ctl文件示例：

------------------------------------------------------------------------------------

dset D:\output\ano.grd  （被描述的文件路径及文件名）

title hight                                        （简略描述数据内容，将在查询命令q中出现）

undef  -9.99e+33                          （说明缺测值）

xdef 144 linear 0 2.500000         （时空维数设定）

ydef 73 linear -90 2.500000

zdef 3 levels 850 500 200

tdef 30 linear 00Z01JUN1979  1mo

vars 1                                              （变量定义）

hgt 0 99 heights

endvars

-------------------------------------------------------------------------------------

说明：

1）时空维数设定

xdef：  x方向格点数360，均匀网格linear， 0（起始经度，负数表示西经），格距2.500000

ydef ： y方向格点数73 ，线性映射linear， -90 （起始经度，负数表示南纬），格距2.500000

zdef： z方向格点数 1 ，levels（顺序枚举出全部等压面）： 850hPa,500hPa，200hPa

tdef ：时次数：30，均匀 linear ，起始时间：00Z01JUN1979 ，时间增量： 1mo

附注：起始时间格式：mm（两位数的分钟）Zdd（一位或两位数的日期）mmm（三字符的月份缩写）yyyy（四位数的年份）；时间增量格式：vv（一位或两位整数）kk（增量类型：分钟mn，小时hr，天dy，mo月，yr年）

 

2）变量定义

var：变量数1

hgt（变量名）， 3（所描述的层次数(易错)，0表示该变量只有一层，且不对应于垂直层，如地表变量）， 99 （若是自己写的ctl写“99“，其他地方的照样摘抄即可）， heights（变量描述，可有可无）

 

Grads 画风场、高度场合成gs

*打开uv风场（1号、2号文件）

'open D:\uc.ctl'

'open D:\vc.ctl'

*打开高度场（3号文件）

'open D:\zc.ctl'

'enable print D:\ForcastNoTimeSmoothing.gmf'

*画等高线(用黑色）

'set ccolor 1'

'd zc.3'

*画uv合成风矢量

'd uc.1;vc.2'

'draw title Forcast'

'draw xlab LON'

'draw ylab LAT'

'print'

'disable print'

;

 

Grads叠加uv画风场

关于画风场: 
先把两个文件都打开比如: 
'open d:/u.ctl' 
'open d:/v.ctl' 
Grads会自动给文件编上号(1，2，3....) 
然后画失量风场要同时display1号文件里的u和2号文件里的v两个变量，所以写成 
'd u.1;v.2'

 

Grads色标问题

增强的色标，实现GrADS颜色标尺间隔标值 ：

http://bbs.06climate.com/forum.php?mod=viewthread&tid=3228

增强色标又多了一种形式-GrADS的不带三角的cbarn色标：

http://bbs.06climate.com/forum.php?mod=viewthread&tid=4587&extra=&page=1

GrADS中自定义色表的使用

http://bbs.06climate.com/forum.php?mod=viewthread&tid=3208&extra=page%3D1

 

Grads 循环出图一例

'sdfopen E:\sst.mnmean.nc'

'set gxout fwrite'

'set fwrite E:\sst-jja.grd'

'set x  1 180'

'set y  1 73'

a=1967

while(a<=2014)

'd ave(sst,time=jun'a',time=aug'a')'

a=a+1

endwhile

'reinit'

 

Grads底图存放及调用

（1）底图存放路径：

根据你安装GRADS的磁盘和路径)\GRADS\Classic\OpenGrADS\Contents\Resources\SupportData\(然后复制中国底图什么的放到这里来)

具体根据你GRADS的版本名字,如果是GRADS2.0.1的版本则变成是：E:(根据你安装GRADS的磁盘和路径)\GRADS.2.0.1\Contents\Resources\SupportData\(然后复制中国底图什么的放到这里来)

拷贝到C:\GrADS19\dat文件夹下

（2）使用方法：

'set mpdset cnworld'

'set mpdset cnscs'

（3）我目前有的底图：cn，cn&river，cnbase，cnriver，cnscs，cnwater，cnworld

 

Grads Skip语句

skip()

 

skip (expr, skipi, skipj)

 

Sets alternating values of expr to the missing data value.

 

expr     A valid grid expression that may have 1 or 2 varying dimensions

skipi    Skip factor in the I dimension of expr

skipj    Skip factor in the J dimension of expr

Usage Notes

 

This function is often used while displaying wind arrows or barbs to thin the number of arrows or barbs. It is not necessary to use the skip function on both the U and V wind components; it is sufficient to populate only one component with missing data values to suppress the plotting of the wind arrow or barb.

Examples

 

Suppose you have a time series of 3-hourly data, but you want to display values at 6-hourly time steps:

 

set x 1

set y 1

set t 1 last

d skip(var,2)

 

To display every other grid point in both the X and Y direction:

 

d skip(u,2,2);v

 

To display every grid point in the Y direction, but every 5th grid point in the X direction:

 

d skip(u,5,1);v

 

This example script "d_uv.gs" written by Wesley Ebisuzaki automatically sets the skip factor based on the plot dimensions. 

*

* This function does a d skip(ugrd,n);v

* where n is automatically set to an appropriate value

*

* usage: d_uv ugrd vgrd 

*

* v1.1 w. ebisuzaki

* v1.2 4/6/98 revised empirical formula for skip

*

function duv(arg)

u = subwrd(arg,1)

v = subwrd(arg,2)

 

* get lat/lon info

'query dims'

lons = sublin(result,2)

lats = sublin(result,3)

dx = subwrd(lons,13) - subwrd(lons,11)

dy = subwrd(lats,13) - subwrd(lats,11)

 

* Determine skip factor 

dn = dx

if (dy > dx) ; dn = dy ; endif

skip = dn / 50 + 0.5

if (skip < 1) ; skip=1 ; endif

 

* Display the plot

'd skip('u','skip');'v