利用站点资料画出等值线图的GrADS示例

由于很多同志是做业务的，接触的都是站点资料，这里将详细的介绍如何利用站点资料，画出某变量的等值线图。

第一步：站点资料（txt格式，即十进制格式）转成二进制格式

GrADS画等值线图是有要求的，即只能画格点上的二进制数据。

首先，将观测资料按如下格式放在txt文件中：

第一列为某站点的纬度，第二列为该站点的经度，第三列为该站点某月平均某一个指数，也就是变量，它可以是降水、气温等。

然后，将该数据转换为二进制，使用以下fortran程序：

c----------------------------------------------------------

c This program is used to change the format of ascii data

c into the form of binary ,which is supported by the GrADS

c----------------------------------------------------------

program main

real, dimension(160) :: lat, lon, ind

open(1,file='d:\\drought\\index.txt',status='old')

do i=1,160

read(1,*) lat(i), lon(i),ind(i)

enddo

close(1)

call stntogrd(ind)

end

subroutine stntogrd(x)

character*8 stid(160)

do i=1,160

stid(i)=char(i)

enddo

open(3,file='e:\\drought\\index.grd',form='binary')

tim=0.0

nlev=1

nflag=1

do i=1,160

write(3) stid(i),lat(i),lon(i),tim,nlev,nflag,x(i)

enddo

nlev=0

write(3) stid(i-1),lat(i-1),lon(i-1),tim,nlev,nflag

close(3)

return

end

于是得到了二进制格式的变量index.grd,接着我们需要给该二进制数据写一个描述性文件如下（新建一个写字板，取名为\"index.ctl\"）：

dset d:/drought/index.grd

dtype station

stnmap d:/drought/index.map

undef -999.0

title drought index

tdef 1 linear Jul1951 1mo

vars 1

ind 0 99 drought index

endvars

到此，二进制数据文件的描述文件写完了，然后，在GrADS命令窗口输入命令“！stnmap” 如：

ga_> !stnmap

系统会提示如：“Enter stn ctl filename:” 输入: d:/drought/index.ctl 系统会生成\"index.map\"文件。第一步完成。

第二步： 制作网格

为什么要做网格呢，因为我们一开始说了，GrADS是画网格点上的二进制数据的，我们现在有了二进制数据，还有网格没有画。利用如下fortran程序制作网格文件(该文件同样也是二进制的)：

Program main

parameter(nx=71, ny=41)

real lat(ny), lon(nx)

real s(nx,ny)

open(1, file= 'd:/drought/grid.grd', form='binary')

lat(1)=15.0

lon(1)=70.0

do j=1,ny-1

lat(j+1)=lat(j)+1.0

enddo

do i=1,nx-1

lon(i+1)=lon(i)+1.0

enddo

do i=1,nx

do j=1,ny

s(i,j)= 1

enddo

enddo

write(1) s

end

运行程序会生成grid.grd的格点文件，接下来为该格点文件写描述性文件，类似于变量的描述性文件：

dset d:/drought/grid.grd

undef -999.0

title Grid data

xdef 71 linear 70 1

ydef 41 linear 15 1

zdef 1 linear 1000 1

tdef 1 linear jul1951 1mo

vars

g 0 99 grid data

endvars

第三步：写脚本

新建一个写字板，取名为 index.gs, 内容可以如下：

'open d:/drough/index.ctl'

'open d:/drought/grid.grd'

'enable print d:/drought/index.gmf'

'set lon 73 135'

'set lat 15 55'

'set mpdset hires cnworld'

'define a=oacres(g, ind.2, 1.5)'

'set gxout shaded'

'set gxout contour'

'd a'

'print'

'disable print'

'reinit'

;

到此结束，在GrADS里面数据命令： run d:/drought/index.gs 就可以画出等值线图了~