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

最新推荐文章于 2021-12-02 16:08:12 发布

hdl62402

最新推荐文章于 2021-12-02 16:08:12 发布

阅读量4.4k

点赞数 1

分类专栏：气象资料

气象资料专栏收录该内容

1 篇文章 0 订阅

订阅专栏

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

由于很多同志是做业务的，接触的都是站点资料，这里将详细的介绍如何利用站点资料，画出某变量的等值线图。
第一步：站点资料（txt格式，即十进制格式）转成二进制格式
         GrADS画等值线图是有要求的，即只能画格点上的二进制数据。
         首先，将观测资料按如下格式放在txt文件中：

https://i-blog.csdnimg.cn/blog_migrate/3a841942d807cede74d407a376399ddc.png (图片地址)
         第一列为某站点的纬度，第二列为该站点的经度，第三列为该站点某月平均某一个指数，也就是变量，它可以是降水、气温等。
         然后，将该数据转换为二进制，使用以下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 就可以画出等值线图了~