GrADS绘图学习技巧与实例(阿木)

以下技巧总结都是笔者从学习实践过程中总结出来的，基本的问题。

不求全面，希望对读者学习有用，如果有问题，敬请留言指正，以促进交流学习！（笔者：阿木）
1、软件综述：grads软件是一款绘图软件除了绘制图形，还可以提取数据，主
要应用是在大气科学中，当然只要是数据处理成grads能够读取的数据文件就可以进行相关绘图。

软件版本问题，软件本身不是很大，我接触到1.8、1.9、
2.0版本的，1.8版本的安装很多情况还要修改环境变量、1.9版本的不识别
‘sdfopen’命令，最稳定的版本是2.0版本，所以笔者推荐学习者安装2.0版本，选择默认安装路径就可以。

2、文件类型简述：grads处理的是网格数据，可以处理的数据类型有：grd、
grib、nc（海洋常用的数据），cdf（雷达卫星数据），其中nc、cdf数据都是自带描述文件，不需要ctl，grib数据要通过命令生成ctl、index数据才可以调用，常用的是grd数据，需要ctl。

3、数据文件转换：grads软件识别的数据是二进制无格式数据，文件类型是
‘binary’，写入和生成时是不需要格式的如read(20) sst(i,j,iz,it),20为文件号，通常是十进制数据与grd数据间转换，这里给一个grd转换成txt数据的fortran程序：
parameter(nx=56,ny=41,nz=1,nt=360)
dimension sst(nx,ny,nz,nt)
real sst
open(15,file='sst.grd',form='binary') ！固定的用form=‘binary’就是二进制数据open(16,file='sst.txt') ！新建txt文件
do it=1,nt
do iz=1,nz
read(15) ((sst(i,j,iz,it),i=1,nx),j=1,ny) ！read后只有文件号，数据是无格式的
enddo
enddo
do it=1,nt
do iz=1,nz
write(16,*) ((sst(i,j,iz,it),i=1,nx),j=1,ny) ！输出时是txt文件可直接看的数据，有格式输出，有*
enddo
enddo
close(15)
close(16)
end
写程序时：注意格点数要与数据对应，如：上程序对应的数据是经度90~200，纬度-20~60，时间：1971.01~2000.12共360个月的海面温度数据，数据格点精度2*2 ，nx=(200-90)/2+1,ny=(60-(-20))/2+1,nt=360，nz=1，大气的数据要根据数据的层次确定几层。

4、grd 、ctl、、gs、nc详述
grd文件：grd数据不可直接看，为二进制无格式数据，简单的说只有1和0，而且数据间没有间隔，grads识别grd是根据ctl进行划分的，根据ctl中的经度、纬度、层次、时间，精度进行数据分块。

ctl实例：
dset C:\data\sst.grd
undef -9.99E+33
title sea surface tempture
xdef 56 linear 90 2
ydef 41 linear -20 2
tdef 360 linear jan1971 1mo
zdef 1 levels 0
vars 1
sst 0 99 surface sea tempture
endvars
上例数据的数据顺序是（以下是数据对应的经纬度）
纬度经度（t=1）
-20 90~200的纬度20S的从90E~200的56个数据
-18 90~200的纬度18S的从90E~200的56个数据
.
.
.
6090~200的纬度60N的从90E~200的56个数据
以上为一层的数据，接下来是t=2，t=3……t=360的数据，每个时间点的每一层是如上格式，
编程时读取和写如的数据循环顺序依次是：时间、层次、纬度、经度，读者参照3中的fortran程序加以理解。

ctl文件：具体其他指导书上都有，我这里强调的是sst 后面的0表示一层，如果是两层以上则是2,3……，1层是固定用法，sst后面的99是默认设置；undef -9.99E+33 此处的数值决定了软件将文件中的那些值认定为不绘制的数据，所以这个值一定要与数据对应。

gs 文件，批命令文件，与ctl一样是用记事本编写，另存为.gs文件，文件都是
命令，方便大段的命令编写、修改，很常用，尤其时绘制的图要求比较多时必须用，免得在命令窗口重复输入命令浪费时间，如下例子：
'reinit'
'open c:\data\uv.ctl'
'set t 7'
'set lev 850'
'set grid off'
'set vpage 0 8.5 0 3.6'
'set parea 0.3 8.5 0.3 3.5'
'set gxout vector'
'set grads off'
'd u;v'
'set vpage 0 8.5 3.6 7.2'
'set grads off'
'set gxout barb'
'd u;v'
'set vpage 0 8.5 7.2 11'
'set grads off'
'set gxout grid'
'd u;v'
'printim c:\images\gxout3uv850.png white'
;
Gs文件以‘；’结尾，分号后面不能有空格，除了循环命令外都需要单引号将命令引起来。

nc文件：此类数据不需要ctl可直接用sdfopen命令直接打开，往往需要知道数据文件中的各个纬度特征，可以用：q ctlinfo命令查询该文件的内置ctl，这样一切都会很清楚，尤其是数据精度：截图如下（下例即2*2的网格经度）：
以上为基础知识介绍，以下为技巧命令
1、绘图时，图的时间下标可用set grads off 命令关掉，网格用set grid off 命令
关掉，需注意的是，网格关一次就一直有效，而下标则是每次绘完图就自动开启，所以建议读者在每次的绘图命令之前加set grads off 。

2、
3、reinit命令是让窗口恢复到刚打开时的界面，会恢复所有的set，同时关掉所
有的打开文件，为为防止前面打开文件的干扰，建议读者在每个gs文件的第一条就加上这条命令（参考gs描述的例子）。

4、
5、数据维度设定：
set lon 90 set x 1
set lat 80 set y 20
set lev 1000 set z 1
set t 1 12 set time jan1970 dec1970
上面两种设定等价，左边是实际维度设定法，右边是给点设定法，读者需要确定具体的格点数。

需要注意的是，在设定全球尺度时，经度0和360是同一个格点，所以set 0 360会出错，这时设定格点的方法比较好：set x 1 180。

6、vpage和parea的区别：
vapge是对整个绘图区分块，需要几张图就划分成几块，给英寸时不需要给标注、标题预留空间，如：
要横着绘制两张图那么就是set vpage 0 5.5 0 8.5
（11*8.5）绘制第一张图
set vpage 5.5 11 0 8.5
绘制第二张图
parea 是描述的绘图时图形的四根边线的大小，而且是虚页的尺寸
如：画一张图set parea 1.0 10.2 0.8 7.8 四个值的范围取决于vapge的长度0<1.0<10.2<11，如果是上例中的两张图，那么：set parea 0.5 4.8 0.8 7.9，0<0.5<4.8<5.5 ，数值范围取决于水平、垂直的长度，与起点无关
如set vpage 0 11 3.5 7
set parea 0.8 10 0.3 3.2
y上满足0<0.3<3.2<(7-3.5) 即可，实际y长度决定范围，与起点无关。

5、'set annot 5 8'
'set xlopts 3 5 0.18'
'set ylopts 3 5 0.18'
上面这两条命令是设定的x、y轴下标数字的颜色、粗细、字号
'set annot 5 8'是设定坐标轴线、标题的颜色粗细，会重置xlopts的部分设定，希望读者注意，命令间的互相干扰
6、cbar cbarn命令
cabr
cbarn
这里只讲解cbarn 命令，
cbarn sf vert xmid ymid
其中sf 为标尺,1为全尺寸0.5为半尺寸；vert为放置位置，0为水平，1为垂直；xmind 、ymid为色标的中心位置
如：cbarn 0.5 0 2.5 1.5 色标半尺度长，水平放置，色标的中心英寸坐标（2.5,1.5）。

如果只是cabrn或cbar命令，会按照默认的全长，水平或数值取决于图的哪边空位大绘制。

强调一点，加色标和加标题都是在绘制出图形后才能加，因为色标是根据阴影图确定对应色值，标题根据图形大小确定标题位置。

7、cmin、cmax是命令是用于绘制大于或小于某数值的线或区域图形，在每次
绘图之后会重设，所以如果多次使用一定要每次d之前加上，这点与set grads off相似。

8、求12个月每个月的海温距平值
'set t 1 12'
'asst=ave(sst,t+1404,t=1764,12)'
'modify asst seasonal'
'set time jan1971'
'd asst'
'set t 1405 1764'
'nasst=sst-asst'
'set time JAN1998'
'd nasst'
set t 1 12 是设定asst变量有12个时间序列，每个时间格点放一个平面的平均值，
asst=ave(sst,t+1404,t=1764,12) 定义变量asst放每个月的平均值，随着t从1变化到12，一次求的每个月的平均值，起始时间是t=1405，终止时间是1764，t=1时，相当于1+13+25+37……，即每年的1月份的值求平均（这里t=1+1404=1405是1971年1月，因此是求1971~2000年每年1月份的平均值），t=2,3.4,，……12 与1同理。

如此将12个月的平均值都放在了asst里。

Modify asst seasonal 本来12个月的平均值知识放在了1~12的时间序列里，（这里调用的文件数据开始时间是1854年1月），所以asst的值只是在设定1~12以内才能画出，假了这条命令，可以将asst的时间序列扩展到所有时间里，使得每一年的每个月对应都是该月30年的平均值，方便后面求距平。

需要注意的，如果set t 1 4 ，那么这条命令的作用是每年的asst的1~4月值是一样的，是该月的平均值。

由于grads本身软件有一定的问题，当你不是设定一段时间而是一个时间点如：set t 1，步、不用modify命令，所有的时间序列都会有该值'set t 1405 1764'
'nasst=sst-asst'
这两句是求30年（1971 01~2000 12）每个月的海温距平值，共360各月的距平值。

这个不难理解，不做解释。

以下为实例：
1、利用所提供的数据文件，绘制出2003年7月60-150E、0-40N区域内700hPa 流线图，且地形高度场超过2000米以上用黑色阴影显示，并且给出相应标题（请包含姓名拼音与学号），最终将图形保存。

所有命令编写于.gs文件中。

'reinit'
'open c:\data\dxgd.ctl'
'open c:\data\uv.ctl'
'set grid off'
'set grads off'
'set map 3 1 5'
'set xlopts 3 5 0.18'
'set ylopts 3 5 0.18'
'set parea 0.7 10.4 0.6 7.8'
'set t 1'
'set lon 60 150'
'set lat 0 40'
'set gxout shaded'
'set rbcols 1 1 1 1 1 1 1'
'set cmin 2000'
'd h'
'set gxout stream'
'set time jul2003'
'set grads off'
'set lev 700'
'd u.2;v.2'
'set annot 5 8'
'draw title DingXiaoli 20081331001'
'printim c:\images\uv200307850.png white'
2、利用所提供的数据文件，画出2002年1-12月120 E、0-40N 200hPa纬向风的纬度-时间剖面图。

图型要求：
（1）纬向风为西风时填色，东风绘制等值线，给出色标，0值线加粗；
（2）X轴标注为“time”,Y轴标注为“lat”,标题标注为“u 100-120E”。

（3）最终将图形保存为gmf格式。

所有命令编写于.gs文件中
'reinit'
'open c:\data\uv.ctl'
'set grid off'
'set grads off'
'set xlopts 3 5 0.16'
'set ylopts 3 5 0.16'
'set parea 1 10.4 1.2 7.9'
'set lon 120'
'set lat 0 40'
'set lev 200'
'set t 1 12'
'set annot 5 8'
'enable print c:\images\uv200.gmf' 'set xyrev on'
'set gxout shaded'
'set cmin 0'
'd u'
'cbarn 1 0 5.5 0.3'
'set grads off'
'set gxout contour'
'set cmax 0'
'set grads off'
'd u'
'set clevs 0'
'set ccolor 2'
'set cthick 10'
'set grads off'
'd u'
'draw xlab time'
'draw ylab lat'
'draw title u 100-120E'
'print'
'disable print'
;
3、十二个月的nc数据整合到一个文件中（fortran）这个程序可用于整合数据parameter(nx=360,ny=181,nz=26,nt=12)
dimension temp(nx,ny,nz,nt)
open(1,file='argo_200901.grd',form='binary')
open(2,file='argo_200902.grd',form='binary')
open(3,file='argo_200903.grd',form='binary')
open(4,file='argo_200904.grd',form='binary')
open(5,file='argo_200905.grd',form='binary')
open(6,file='argo_200906.grd',form='binary')
open(7,file='argo_200907.grd',form='binary')
open(8,file='argo_200908.grd',form='binary')
open(9,file='argo_200909.grd',form='binary')
open(10,file='argo_200910.grd',form='binary')
open(11,file='argo_200911.grd',form='binary')
open(12,file='argo_200912.grd',form='binary')
open(37,file='argo_2009.grd',form='binary')
do iz=1,nz
read(1) ((temp(i,j,iz,1),i=1,nx),j=1,ny)
read(2) ((temp(i,j,iz,2),i=1,nx),j=1,ny)
read(3) ((temp(i,j,iz,3),i=1,nx),j=1,ny)
read(4) ((temp(i,j,iz,4),i=1,nx),j=1,ny)
read(5) ((temp(i,j,iz,5),i=1,nx),j=1,ny)
read(6) ((temp(i,j,iz,6),i=1,nx),j=1,ny)
read(7) ((temp(i,j,iz,7),i=1,nx),j=1,ny)
read(8) ((temp(i,j,iz,8),i=1,nx),j=1,ny)
read(9) ((temp(i,j,iz,9),i=1,nx),j=1,ny)
read(10) ((temp(i,j,iz,10),i=1,nx),j=1,ny)
read(11) ((temp(i,j,iz,11),i=1,nx),j=1,ny)
read(12) ((temp(i,j,iz,12),i=1,nx),j=1,ny)
enddo
do it=1,nt
do iz=1,nz
write(37) ((temp(i,j,iz,it),i=1,nx),j=1,ny)
enddo
enddo
close(1)
close(2)
close(3)
close(4)
close(5)
close(6)
close(7)
close(37)
end
4、竖着绘制三张不同输出形式的uv风场图（parea可保证图形达到尽可能的大）'reinit'
'open c:\data\uv.ctl'
'set t 7'
'set lev 850'
'set grid off'
'set vpage 0 8.5 0 3.6'
'set parea 0.3 8.5 0.3 3.5'
'set gxout vector'
'set grads off'
'd u;v'
'set vpage 0 8.5 3.6 7.2'
'set grads off'
'set gxout barb'
'd u;v'
'set vpage 0 8.5 7.2 11'
'set grads off'
'set gxout grid'
'd u;v'
'printim c:\images\gxout3uv850.png white'
;
5、用循环绘制矩形块，是两个起点逆时针绘制的命令，这个不需要数据文件。

'reinit'
k1=-0.8
m1=-0.3
j=1
while(j<=10)
j=j+1
k1=k1+1.2
m1=m1+1.2
k2=10.4-k1
m2=11.4-m1
'set line 'j' 2 6'
'draw recf 'k1' 0.2 'm1' 0.7'
'draw recf 'k2' 7.8 'm2' 8.3'
endwhile
k1=-0.75
m1=-0.25
j=1
while(j<=10)
j=j+1
k1=k1+0.95
m1=m1+0.95
k2=8.0-k1
m2=9.0-m1
'set line 'j' 2 6'
'draw recf 10 'k1' 10.5 'm1''
'draw recf 0.4 'k2' 0.9 'm2''
endwhile
;
11。