GrADS绘图学习技巧与实例.docx
《GrADS绘图学习技巧与实例.docx》由会员分享，可在线阅读，更多相关《GrADS绘图学习技巧与实例.docx（14页珍藏版）》请在冰点文库上搜索。
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）
dimensionsst（nx,ny,nz,nt）
realsst
open（15,file='sst.grd',form='binary'）！
固定的用form=‘binary’就是二进制数据
open（16,file='sst.txt'）！
新建txt文件
doit=1,nt
doiz=1,nz
read（15）（（sst（i,j,iz,it）,i=1,nx）,j=1,ny）！
read后只有文件号，数据是无格式的
enddo
enddo
doit=1,nt
doiz=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实例：
dsetC:
\data\sst.grd
undef-9.99E+33
titleseasurfacetempture
xdef56linear902
ydef41linear-202
tdef360linearjan19711mo
zdef1levels0
vars1
sst099surfaceseatempture
endvars
上例数据的数据顺序是（以下是数据对应的经纬度）
纬度经度（t=1）
-2090~200的纬度20S的从90E~200的56个数据
-1890~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'
'openc:
\data\uv.ctl'
'sett7'
'setlev850'
'setgridoff'
'setvpage08.503.6'
'setparea0.38.50.33.5'
'setgxoutvector'
'setgradsoff'
'du;v'
'setvpage08.53.67.2'
'setgradsoff'
'setgxoutbarb'
'du;v'
'setvpage08.57.211'
'setgradsoff'
'setgxoutgrid'
'du;v'
'printimc:
\images\gxout3uv850.pngwhite'
;
Gs文件以‘；’结尾，分号后面不能有空格，除了循环命令外都需要单引号将命令引起来。
nc文件：
此类数据不需要ctl可直接用sdfopen命令直接打开，往往需要知道数据文件中的各个纬度特征，可以用：
qctlinfo命令查询该文件的内置ctl，这样一切都会很清楚，尤其是数据精度：
截图如下（下例即2*2的网格经度）：
以上为基础知识介绍，以下为技巧命令
1、绘图时，图的时间下标可用setgradsoff命令关掉，网格用setgridoff命令关掉，需注意的是，网格关一次就一直有效，而下标则是每次绘完图就自动开启，所以建议读者在每次的绘图命令之前加setgradsoff。
2、
3、reinit命令是让窗口恢复到刚打开时的界面，会恢复所有的set，同时关掉所有的打开文件，为为防止前面打开文件的干扰，建议读者在每个gs文件的第一条就加上这条命令（参考gs描述的例子）。
4、
5、数据维度设定：
setlon90setx1
setlat80sety20
setlev1000setz1
sett112settimejan1970dec1970
上面两种设定等价，左边是实际维度设定法，右边是给点设定法，读者需要确定具体的格点数。
需要注意的是，在设定全球尺度时，经度0和360是同一个格点，所以set0360会出错，这时设定格点的方法比较好：
setx1180。
6、vpage和parea的区别：
vapge是对整个绘图区分块，需要几张图就划分成几块，给英寸时不需要给标注、标题预留空间，如：
要横着绘制两张图那么就是setvpage05.508.5
（11*8.5）绘制第一张图
setvpage5.51108.5
绘制第二张图
parea是描述的绘图时图形的四根边线的大小，而且是虚页的尺寸
如：
画一张图setparea1.010.20.87.8四个值的范围取决于vapge的长度0<1.0<10.2<11，如果是上例中的两张图，那么：
setparea0.54.80.87.9，
0<0.5<4.8<5.5，数值范围取决于水平、垂直的长度，与起点无关
如setvpage0113.57
setparea0.8100.33.2
y上满足0<0.3<3.2<（7-3.5）即可，实际y长度决定范围，与起点无关。
5、'setannot58'
'setxlopts350.18'
'setylopts350.18'
上面这两条命令是设定的x、y轴下标数字的颜色、粗细、字号
'setannot58'是设定坐标轴线、标题的颜色粗细，会重置xlopts的部分设定，希望读者注意，命令间的互相干扰
6、cbarcbarn命令
cabr
cbarn
这里只讲解cbarn命令，
cbarnsfvertxmidymid
其中sf为标尺,1为全尺寸0.5为半尺寸；vert为放置位置，0为水平，1为垂直；xmind、ymid为色标的中心位置
如：
cbarn0.502.51.5色标半尺度长，水平放置，色标的中心英寸坐标（2.5,1.5）。
如果只是cabrn或cbar命令，会按照默认的全长，水平或数值取决于图的哪边空位大绘制。
强调一点，加色标和加标题都是在绘制出图形后才能加，因为色标是根据阴影图确定对应色值，标题根据图形大小确定标题位置。
7、cmin、cmax是命令是用于绘制大于或小于某数值的线或区域图形，在每次绘图之后会重设，所以如果多次使用一定要每次d之前加上，这点与setgradsoff相似。
8、求12个月每个月的海温距平值
'sett112'
'asst=ave（sst,t+1404,t=1764,12）'
'modifyasstseasonal'
'settimejan1971'
'dasst'
'sett14051764'
'nasst=sst-asst'
'settimeJAN1998'
'dnasst'
sett112是设定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里。
Modifyasstseasonal本来12个月的平均值知识放在了1~12的时间序列里，（这里调用的文件数据开始时间是1854年1月），所以asst的值只是在设定1~12以内才能画出，假了这条命令，可以将asst的时间序列扩展到所有时间里，使得每一年的每个月对应都是该月30年的平均值，方便后面求距平。
需要注意的，如果sett14，那么这条命令的作用是每年的asst的1~4月值是一样的，是该月的平均值。
由于grads本身软件有一定的问题，当你不是设定一段时间而是一个时间点如：
sett1，步、不用modify命令，所有的时间序列都会有该值
'sett14051764'
'nasst=sst-asst'
这两句是求30年（197101~200012）每个月的海温距平值，共360各月的距平值。
这个不难理解，不做解释。
以下为实例：
1、利用所提供的数据文件，绘制出2003年7月60-150E、0-40N区域内700hPa流线图，且地形高度场超过2000米以上用黑色阴影显示，并且给出相应标题（请包含姓名拼音与学号），最终将图形保存。
所有命令编写于.gs文件中。
'reinit'
'openc:
\data\dxgd.ctl'
'openc:
\data\uv.ctl'
'setgridoff'
'setgradsoff'
'setmap315'
'setxlopts350.18'
'setylopts350.18'
'setparea0.710.40.67.8'
'sett1'
'setlon60150'
'setlat040'
'setgxoutshaded'
'setrbcols1111111'
'setcmin2000'
'dh'
'setgxoutstream'
'settimejul2003'
'setgradsoff'
'setlev700'
'du.2;v.2'
'setannot58'
'drawtitleDingXiaoli20081331001'
'printimc:
\images\uv200307850.pngwhite'
2、利用所提供的数据文件，画出2002年1-12月120E、0-40N200hPa纬向风的纬度-时间剖面图。
图型要求：
（1）纬向风为西风时填色，东风绘制等值线，给出色标，0值线加粗；
（2）X轴标注为“time”,Y轴标注为“lat”,标题标注为“u100-120E”。
（3）最终将图形保存为gmf格式。
所有命令编写于.gs文件中
'reinit'
'openc:
\data\uv.ctl'
'setgridoff'
'setgradsoff'
'setxlopts350.16'
'setylopts350.16'
'setparea110.41.27.9'
'setlon120'
'setlat040'
'setlev200'
'sett112'
'setannot58'
'enableprintc:
\images\uv200.gmf'
'setxyrevon'
'setgxoutshaded'
'setcmin0'
'du'
'cbarn105.50.3'
'setgradsoff'
'setgxoutcontour'
'setcmax0'
'setgradsoff'
'du'
'setclevs0'
'setccolor2'
'setcthick10'
'setgradsoff'
'du'
'drawxlabtime'
'drawylablat'
'drawtitleu100-120E'
'print'
'disableprint'
;
3、十二个月的nc数据整合到一个文件中（fortran）这个程序可用于整合数据
parameter（nx=360,ny=181,nz=26,nt=12）
dimensiontemp（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'）
doiz=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
doit=1,nt
doiz=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'
'openc:
\data\uv.ctl'
'sett7'
'setlev850'
'setgridoff'
'setvpage08.503.6'
'setparea0.38.50.33.5'
'setgxoutvector'
'setgradsoff'
'du;v'
'setvpage08.53.67.2'
'setgradsoff'
'setgxoutbarb'
'du;v'
'setvpage08.57.211'
'setgradsoff'
'setgxoutgrid'
'du;v'
'printimc:
\images\gxout3uv850.pngwhite'
;
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
'setline'j'26'
'drawrecf'k1'0.2'm1'0.7'
'drawrecf'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
'setline'j'26'
'drawrecf10'k1'10.5'm1''
'drawrecf0.4'k2'0.9'm2''
endwhile
;
升级会员 