用Grads画等值线（二）------.map生成

      要生成.map文件，需要以下四步： （一）准备一个.ch文件，其中有需要的站点数据； （二）编写一个简单的程序把.ch中的数据写成二进制格式的.dat文件。（三）编写一个站点控制文件.ctl。（四）运行如下Grads命令：stnmap -i C:/Grads/rain.ctl，其中"C:/Grads/rain.ctl"指站点控制文件的路径及名称。

     上面生成.map文件的例子在论坛中也有不少，大都是用fortran写的程序，也有少数用C写的。但是我自己写的到了第四步总是报错。我相信对刚学Grads而且非气象专业的人来说，也是常碰到的。可是对熟悉Grads的人来说可能这根本不算个问题吧，所有很少有人说，而且说到的说得也不够详细透彻。后来费了九牛二虎之力终于弄明白了。站点数据有非常严格的格式，格式是按如下规定的：

 

    　站点资料采用二进制文件进行存储, 按每个时次一个报告的顺序直接记录各站资料, 每个时次的站点报告组又分成两部分: 地面报和高空报。站点报数据的排放顺序为:
　　① 提供站点经纬度的头记录
　　② 地面报变量
　　③ 高空报变量
　　④ 一个时次完成后加上一个特殊头记录 (没有数据组) , 表示本时次数据报结束。
      头记录共占28 个字节, 包括:
     STA T: 站点标识, 字符型 (st r ing 或char) , 占8个字节; 其赋值可以是任意数字, 也可以是字母。
     RLON: 经度, 单精度型 (single) , 占4 个字节。
     RLA T: 纬度, 单精度型 (single) , 占4 个字节。
     T I M : 本报的时次, 单精度型(single) , 占4 个字节,
     NL EV: 头记录后的数据组数, 单精度型 (single) ,占 4 个字节, 包括一个地面报组加上各高度层数, 取0 时标志一个时次报文的结束。
     N FLA G: 单精度型 (single) , 占4 个字节。取0 表示头记录的报文中没有地面变量, 取1 时则有。
　　在头记录之后, 顺序记录本时次的各报文, 首先写入地面变量, 然后是高空变量。当本时次的所有报告写完后, 写出一个特别记录, 该记录没有数据组数, 以表明一个时间组的结束, 然后是下一时次记录。

      不仅要求数据顺序，而且每个字段所占字节数都是固定的，所有不知道而任写很容易出错。.ch文件中的数据格式也要用空格和换行符严格控制，所有最好用UIEditor等编辑器来编辑。下面是我所使用的营口气象站的雨量数据：

 L6001  40.674  122.328 144.0
 L6002  40.829  122.175 105.0
 L6016  40.530  122.621  44.0
 L6017  40.568  122.684  63.0
 L6019  40.481  122.864  32.0
 L6003  40.850  122.277 129.0
 L6004  40.861  122.434 124.0
 L6005  40.780  122.385 126.0
 L6006  40.690  122.561  96.0
 L6015  40.528  122.462  74.0
 L6020  40.493  122.796  38.0
 L6010  40.634  122.391 123.0
 L6000  40.674  122.273 126.0
 L6007  40.686  122.380 123.0
 L6009  40.596  122.266 138.0
 L6036  40.432  122.610  40.0
 L6011  40.599  122.516  75.0
 L6021  40.397  122.875  49.0
 L6082  40.184  122.018  50.0
 L6053  40.108  122.002  50.0
 L6047  40.139  122.187  51.0
 L6038  40.338  122.393  40.0
 L6042  40.217  122.144  82.0
 L6012  40.577  122.476 144.0
 L6008  40.338  122.526 101.0
 L6032  40.383  122.504  46.0
 L6041  40.256  122.734  68.0
 L6045  40.211  122.519  91.0
 L6046  40174   122.706  27.0
 L6050  40.095  122.109  50.0
 L6030  40.493  122.427  76.0
 L6031  40.416  122.655  82.0
 L6033  40.396  122.724  37.0
 L6034  40.343  122.236  97.0
 L6035  40.433  122.472  37.0
 L6044  40.243  122.439  86.0
 L6048  40.231  122.249  60.0
 L6049  40.134  122.545  94.0
 L6051  40.063  122.493 104.0
 L6052  40.101  122.248  48.0
 L6037  40.263  122.124  83.0
 L6040  40.269  122.186  70.0
 L6083  40.399  122.209 129.0

 

 生成rain.dat文件的fortran代码如下：

 program YKRain parameter (st=43) ! st means nums of station character stid(st)*8 integer nflag,nlev real std(st),tim real lat(st),lon(st) open(30,file='D:/YK/rain.ch') open(40,file='D:/YK/rain.dat',form='binary') do i=1,st read(30,'(a8,f8.3,f8.3,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(43),lat(43),lon(43),tim,nlev,nflag close(40) end

 

     rain.ctl的文件如下：

                   dset D:/YK/rain.dat
                   dtype station
                   stnmap D:/YK/rain.map
                   undef -999.0
                   title rainfall
                   tdef 1 linear jan2010 1mo
                  vars 1
                  srainval 0 999.9 rainfall
                  endvars

     如果生成map文件错误，反复查找，看看是不是数据格式的问题，相信一定可以正确生成。