利用站点资料画出等值线图的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 就可以画出等值线图了~

mofangbao

传说中的谁 发表于 2011-6-27 16:34

                               
登录/注册后可看大图

另外这种方法应该配合maskout或province-basmap一起用，避免超出地理边界的情况出现。

恩，你说的对，mask和basemap对于熟悉的人很简单，对于新手还是比较头晕的，如果您有空的话，欢迎您发一个专门的帖子来介绍下哈，谢谢支持喽

传说中的谁

嘟嘟鱼 发表于 2011-6-27 17:28

                               
登录/注册后可看大图

插值半径不能取的太大，否则即使插出来，也是不真实的值，这个需要你不断的尝试然后再跟实际情况对比，找 ...

谢谢版主大大，看来这个问题只有在用的时候不断的调了。

wxxxren

artmunich 发表于 2011-7-5 22:59

                               
登录/注册后可看大图

由于我只是想上手，就用鱼楼主给的部分数据来做test
数据就是楼主的了，命名为test.txt，下面先些了te ...

斑竹，我也遇到这样的问题了。画出来的图和你的都样，错误信息提示也一样，怎么解决的啊，求解！急急急！不胜感激！

mofangbao

很好，可以作为我们自己的新手教程:lol

/db←大哥

这个好！:lol

zhaolanxiai

我按照嘟嘟鱼的程序，怎么还是弄不出来呢？:Q

传说中的谁

可否说明一下cressman插值方法中，插值半径对最后结果的影响？

传说中的谁

另外这种方法应该配合maskout或province-basmap一起用，避免超出地理边界的情况出现。

mofangbao

zhaolanxiai 发表于 2011-6-27 16:19

                               
登录/注册后可看大图

我按照嘟嘟鱼的程序，怎么还是弄不出来呢？

可以把你的问题详细的发个帖子哈，嘿嘿

mofangbao

传说中的谁 发表于 2011-6-27 16:32

                               
登录/注册后可看大图

可否说明一下cressman插值方法中，插值半径对最后结果的影响？

这个改了试下就知道啦，好像是这样的，比如降水，如果插值半径选的不合适，会出现<0的降水量的

嘟嘟鱼

谢谢管理员的解释，我在这里再说明一下，这个做图过程只是一个基本过程，在真正使用的过程中，不能完全照搬的，要根据自己的具体情况来做，gs文件中的过程也只是最基本的省去了许多对图形的设置，这些需要大家在做具体问题的过程中慢慢摸索，摸索也是一个很有趣的过程！！