用grads画站点资料的等值线
一 先把站点文件转化为二进制文件
假定站点资料为
station  lat  lon  var
….   …. ….   …..   
  1 单个时间变量，单个要素 bintd.f90
program bintd
         parameter (st=3276) ! st means nums of station
         character stid(st)*8
     integer nflag,nlev
         real lat(st),lon(st),std(st),tim
         open(30,file='02070108.txt')
         open(40,file='td.dat',form='binary')
         do i=1,st
           read(30,'(a5,2f8.2,f7.1)')stid(i),lat(i),lon(i),std(i)
           write(*,*)stid(i),lat(i),lon(i),std(i)
         end do 
         close(30)
         tim=0.0
         nflag=1
         nlev=1
         do i=1,st
          write(40)stid(i),lat(i),lon(i),tim,nlev,nflag,std(i)
         end do
         nlev=0
         write(40)stid(3276),lat(3276),lon(3276),tim,nlev,nflag
         close(40)
         end
  2 多个时间变量，单个要素binstat.f90
假定站点资料为
time  station  lat  lon  var
……. ……    …  …   ….
program binstat
         character stid*8
     integer nflag,nlev
         real lat,lon,std,tim
         open(30,file='st020120.txt')
         open(40,file='bist020120.dat',form='binary')
         iflag=0
         i=0
10       read(30,'(i3,1x,a5,2f8.2,f7.1)',end=90)id,stid,lat,lon,var
        
       
         if(iflag.eq.0)then
           iflag=1
           ims=im
           ids=id
         end if
         if(ids.ne.id)then
            nlev=0
            write(40)stid,lat,lon,tim,nlev,nflag
            i=i+1
            write(*,*)i
         end if
         ims=im
         ids=id
         tim=0.0
         nflag=1
         nlev=1
         write(40)stid,lat,lon,tim,nlev,nflag
         write(40)var
         go to 10
90       continue
         nlev=0
         write(40)stid,lat,lon,tim,nlev,nflag
         stop
         end    
二 创建站点文件的描述文件 
  1单个时间变量 单个要素  td.ctl
dset td.dat
dtype station
stnmap td.map
undef 999.9
title daily td
tdef 1 linear Jan2002 1mo
vars 1
td 0 99 tdpoint data
endvars
  2 多个时间变量 单个要素  stat.ctl
dset bist020120.dat
dtype station
stnmap svar020120.map
undef 999.9
title daily svar
tdef 31 linear Jan2002 1mo
vars 1
svar 0 99 svars data
endvars
三 stnmap –i xxx.ctl
四 生成相应的格点文件   wgrid.f90
program wgrid
        ! to get  grids
        real gtd(50,50)
        ! redefine the grids numbers when necessary
        open(10,file='/Models/czh/wgrid.dat',form='binary')
        do i=1,50
          do j=1,50
           gtd(i,j)=0
          end do
        end do
        write(10)((gtd(i,j),i=1,50),j=1,50)
        end
五 生成格点文件的描述文件 grid.ctl
dset wgrid.dat
title gtd
undef -9.99E33
xdef 50 linear 115.07 0.18
ydef 50 linear 39.29 0.18
zdef 1 linear 1 1
tdef  1 linear jan2002 1mo
vars 1
gtd 0 99 grid td
endvars
六 生成最后的gs文件 stat.gs
1 显示一张图
'reinit'
'open gtd.ctl'
'open td.ctl'
'set mpdset cnworld'
'set lon 102.2 125.8'
'set lat 31.4 46.4'
'set t 1'
'define tds=oacres(gtd,std.2)'
'set cint 0.5'
'set gxout shaded'
'd tds'
'set gxout contour'
'd tds'
'printim td.gif white'
;
2 同时显示多个时刻的图
'reinit'
'open gtd.ctl'
'open 200201.ctl'
'set mpdset cnworld'
'set lon 115.17 117.672'
'set lat 39.29 41.216'
tm=1
while(tm<=31)
'set t '%tm
'define statvar=oacres(gtd,svar.2)'
'set cint 0.5'
'set gxout shaded'
'd statvar'
'set gxout contour'
'd statvar'
'printim statvar'% tm % '.gif white'
'd statvar'
'c'
tm=tm+1
endwhile

原文出处：http://bbs.lasg.ac.cn/bbs/archiver/tid-10759.html

！

！＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

！

最近在使用站点资料，
一直不满意grads对站点资料处理的局限
在网上查了好多帖子和资料,发现并不是原来我认为的那么局限!
有一些收获总结一下：
1.用fortran写站点资料
有很多的帖子都在讨论这个问题，相信大家都清楚格式怎么写
我只说一个问题：区站号的写法
好多人都是随便写8个字符作为区站号写入文件，当然如果用不到区站号的话就无所谓了；
但是当想用区站号作为标识的时，这样就不合适了，需要认真的写入如‘57073’这样的区站号。
有帖子中提到这8个字符应该怎么排呢？
是
'   57073'
还是
'57073   '
'12345678'
?
应该选择后者,原因是:虽然区站号需要8个字符读入,但实际上只有前7位是有效的,最后一位是作为分隔符的(这个可以在 用set gxout findstn之后的变量属性看出)
为了以后处理方便,我就这样写了,'57073xx'
character(len=5) sta !读入的区站号
character(len=8) id
write(id,"(a5,'xx',1x)")sta

2.使用站点资料中的单站时间序列
'reinit'
'open *.ctl'
*打开站点资料
'set x 1'
'set y 1'
'set z 1'
'set t 1 120'
*设定一维可变
'set stid on'
'd r(stid="57073xx")'
这时你就会发现r变量在站点'57073xx'的时间序列就画出来了!
===============================================
set stid
set stid on|off
Controls whether the station id is displayed next to the values or not.
以前一直不清楚这个是干什么的,现在明白了它表示
变量属性中 stid这个标识的开关
我是很偶然的想明白的,用站点资料试验varname(stid=ident)时,出现错误提示:
Invalid dimension expression
Expecting X/Y/Z/T/LON/LAT/LEV/TIME, found stid
突然明白这是需要一个特殊标识的申明
结论,varname(a=ident)
a 只可以是X/Y/Z/T/LON/LAT/LEV/TIME/stid 其中的一个.
====================================
现在得到了单站的时间序列,接下来就可以计算好多东西
就可以避免在fortran和grads直接来回输出输入数据的麻烦了
但是还是没有格点数据那么方便,函数使用有限制........
http://bbs.lasg.ac.cn/cgi-bin/forum/topic.cgi?forum=3&topic=6583
引用:
下面引用由funny在 2007/03/22 01:17pm 发表的内容：
2) 站点分析
从格点插到站点gr2stn(grid_expr, { stn_expr | lon, lat }, -a )
收集生成格点场coll2gr(cnum, { num | -u } )
Cressman插值  oacres(grid_expr, stn_expr, radii, 1st_guess)
网格平均      oabin(grid_expr, stn_expr, { -f | -c } )
时间平均      stnave(expr, dexpr1, dexpr2, -m count)
时间序列最小值stnmin(expr, dexpr1, dexpr2, -m count)
时间序列最大值stnmax(expr, dexpr1, dexpr2, -m count)
转成1-D格点序列s2g1d(expr)
所有数学函数  abs,cos,acos,sin,asin,tan,atan2,exp,log,log10,pow,mag,sqrt
看起来s2g1d(expr) 可以解决这个问题,可惜我的版本上没有这个函数?
大家试试了,告诉我结果啊!
引用:
s2g1d
s2g1d(expr)
This function converts a station data time series into a 1-dimensional grid. This allows more graphics and analytical operations.
expr            A valid GrADS 1D station data expression in which time is the only varying dimension.

Usage Notes
Examples
'set gxout bar'
'set barbase 0'
'set bargap 10'
'display s2g1d(precip(stid=kdca))'
'set gxout contour'
'display smth9(s2g1d(t(stid=kord)))'
'set gxout vector'
'display const(s2g1d(u(stid=ksux)),0);s2g1d(u(stid=ksux));s2g1d(v(stid=ksux))'

3. set gxout findstn 的使用
'reinit'
'open d:r160outsta-r.ctl'
'set lon 73 136'
'set lat 17 54'
'd r'
'q gxinfo'
say result
'c'
'set gxout findstn'
'd r;5;4'
*这个x,y是窗口中的坐标位置，前一个命令可以查看坐标范围
*然后显示stid, lon,lat,value
id=subwrd(result,1)
x=subwrd(result,2)
y=subwrd(result,3)
value=subwrd(result,4)
say '========================'
say 'stid: ' id
say 'value: ' value
say 'lon: ' x
say 'lat: ' y
'set lon 'x''
'set lat 'y''
'set z 1'
'set t 1 120'
*必须有一维是要变化的
'set stat on'
'd r'
say '========================'
say result
'set stat off'
==================================================
我想说明的是,
可能产生的对
'set gxout findstn'
'd r;x;y'
输入和输出量的误解
*这个x,y是窗口中的坐标位置
*而屏幕上显示的result是stid, lon,lat,value
====================================================================
The graphics output type findstn requires three arguments to be provided with the display command.
The first argument is a station data expression. The 2nd and 3rd arguments are the X and Y screen coordinates of the of the desired search location.
GrADS will find the station closest to the specified X and Y position, and print its stid, lon, and lat.
This graphics output type should only be used when X and Y are the varying dimensions and AFTER a regular display command (that results in graphics output) is entered.

4 oacres 函数的参数设置
引用:
setting the first guess in oacres
http://atm.ucdavis.edu/program/computing_facilities/grads/gaup151.htm
The built-in GrADS obs->grid analysis function "oacres" sets the initial value
of the analysis grid to the arithmetic average of the obs in the area. For po
sitive definate quantities like precipitation, this can produce an unrealistic
analysis in regions of NO obs (e.g., no rain is better guess than average rai
n).
The call to oacres (stands for Cressman r**2 scan analysis) has the form,
d oacres(grid,obs,r1,r2,r3,...,r30)
where,
grid = is a grid to which the obs will be analyzed to obs =is the obs "grid" r
1,...,r30 are the scan radii in GRID UNITS.
The default is:
r1=10.0
r2=7.0
r3=4.0
r4=2.0
r5=1.0
or five radii. This is good for meteorological fields, but may not yield a des
irable analysis for hydrologic fields which are not as continuous. To change t
he first guess, set the penultimate r to -1 and the last r to the desired firs
t guess. For example,
d oacres(pr.1,pr.2,5,4,-1,-0.01)
would do an analysis of the pr.2 obs to the pr.1 grid with 2 scan radii of 5 a
nd 4 grid units with a first guess of -0.01. oacress does not abort on errors
in the parameters. Rather, the default oacres processing is applied.
引用:
下面引用由沙鸥在 2008/03/05 11:22am 发表的内容：
倒数第二个参数为-1，表示倒数第一个参数为做插值时的初猜值，也就是插值前每个格点上的给它赋予的一个数值。如果不设置初猜场，即后面两个参数“－1”和“－0。01”不设置的话，则插值程序会根据站点资料去算一个平均值作为它的初猜值。
常用的气象要素多为连续场,直接使用
d oacres(grid,obs)
默认表示
d oacres(grid,obs,10.0,7.0,4.0,2.0,1.0)
从插值的格点向外搜索10个,7个...格距半径内的站点,
根据这些搜索范围的站点分布密度计算平均值,并且比较,决定采用那种半径, 
最后插值格点的数值由平均值得出.
当然可以自己设置搜索半径
d oacres(grid,obs,r1,r2,r3,...,r30)
具体算法没有仔细看
见
The Cressman Analysis scheme is described in a 1959 paper in Monthly Weather Review
an operational  objective analysis system. GEORGE P. CRESSMAN( Monthly Weather Review).pdf
------------------------------------------------
还有一些数据是离散分布的,水文和地质上常常遇到,这也是克里格分析法的所长之处
d oacres(pr.1,pr.2,5,4,-1,-0.01)
在搜索半径之后多加了两位,
倒数第二位为'-1',标识要使用最后一位那个数作为要插值格点的初始数据,当然这是格式规定啦,程序没有搜索到这个'-1'就不会进行这个子模块啦
倒数第一位为要插值格点的初始数据,它和搜索之后得到的平均值比较,然后做取舍

原文出处：http://bbs.lasg.ac.cn/bbs/thread-6678-1-8.html