micaps 地面观测资料画6小时降水

luoziwuhui

最近，用micaps 地面观测资料画6小时降水的图。在此过程中，参考了清风大大的帖子（GrADS站点资料作图详细解决方案）：@mofangbao
http://bbs.06climate.com/forum.php?mod=viewthread&tid=4903&extra=page%3D2
以及@平流层的萝卜 大大的帖子（micaps多时次单变量批处理及其后grads画图）：http://bbs.06climate.com/forum.php?mod=viewthread&tid=13422 只不过他的是画能见度，我的是画降水的，都差不多，但还是有一点区别。还是贴出来给大家看看吧。高手可以路过啦！

最后出图的时候，还遇到了与@wode_q_x 相同的问题（绘制非全国地图时，用cnbasemap出现错误 ）：http://bbs.06climate.com/forum.php?mod=viewthread&tid=12537
幸好这些问题都得到了接近。跟大家分享一下：

我用到的数据是Micaps 地面观测资料（surface中的plot文件夹下的数据）。如下：

这是第一类数据格式：用于地面填图，变量分别为：

区站号  经度  纬度  拔海高度  站点级别  总云量  风向  风速  海平面气压（本站气压）  3小时变压  过去天气1  过去天气2  6小时降水  低云状  低云量  低云高  露点  能见度  现在天气  温度 中云状  高云状  船向  船速

降水是第13个变量。

   我用到的时次为2012年7月20日02时到23日23时，共32个时次的数据。

详细过程：

1，首先，把数据文件名写到一个文件夹中，用Windows 的CMD命令。打开运行-输入cmd ，回车，然后cd /d e:\rain\plot(我的数据所在的目录)，然后dir/d *.000>filename.txt把该目录下的所有.000文件名写到一个filename.txt中方便下面读取。（感谢清风大大，http://bbs.06climate.com/forum.php?mod=viewthread&tid=3079）

2，用FORTRAN 读取数据，程序如下；

       program read_surface observation
          parameter(nt=32)
          integer::n,num,n1,n2,n3,n4
           real,allocatable::lon(:),lat(:),rd(:)    !动态数组分别用来储存站点的降水
           real:: temp(9)
          character*8,allocatable::sta(:)      !动态数组，储存站号，站号必须是8个字符，5个的话会stnmap出不出来   
        character*12 filename(nt)        !用于储存32个时次的文件名

           open(1,file='E:\rain\plot\filename.txt')
           do i=1,nt
               read(1,*) filename(i)
               print*,'Filename:',filename(i)
            enddo
            close(1)               
        open(2,file='E:\rain\plot\rain.grd',form='binary')  !储存32个时次的降水
       do k=1,nt
        open(3,file=filename(k))
        read(3,*)
        read(3,*) n1,n2,n2,n4,n    !将该时次的站点数赋值于n
        allocate(lat(n))  ！分配动态数组
        allocate(lon(n))
        allocate(rd(n))
        allocate(sta(n))
        close(3)
        open(4,file=filename(k))   !获知n后，重新读取数据。
        read(4,*)
        read(4,*)
          do i=1,n        !只读降水，别的要素暂略
             read(4,*) sta(i),lon(i),lat(i),(temp(j),j=1,9),rd(i)
              tim=0.0;nlev=1;nflag=1
              write(2)  sta(i),lat(i),lon(i),tim,nlev,nflag,rd(i)  !一个时次的降水场输入完毕
           enddo
           nlev=0
            write(2) sta(n),lat(n),lon(n),tim,nlev,nflag
        !一个时次的降水场输入完毕，告诉Grads 该时次的数据结束。这是一个特殊标记。
        deallocate(lat)
        deallocate(lon)
        deallocate(rd)
        deallocate(sta)

！释放动态数组
        close(4)

        enddo
         close(2)
      end program read_surface observation

3，为该数据编写ctl文件，

dset    E:\rain\plot\rain.grd
dtype   STATION
stnmap  E:\rain\plot\rain.map
undef   9999.
title   RAINFALL ONSET DATE IN QTP
tdef 32 linear 02z20Jul2012 3hr
vars    1
r      0  99  Rainfall date
endvars
4， 生成格点背景场（详见mofangbao大大的帖子，在此不罗嗦 ），我用的是0.5*0.5的格点。

为格点背景场写ctl 数据描述文件的时候注意时次一定与降水的数据文件一致（tdef 32 linear 02z20Jul2012 3hr
）。

5， Grads 画图，程序如下：

'reinit'
'c'
'open E:\rain\plot\prc16.ctl'
'open E:\rain\plot\201207r3.ctl'

'enable print E:\rain\plot\201207r3.gmf'

**'set lat 20 40'
**'set lon 85 110'
**'set t 1 32'
**'define rr=oacres(prc,r.2)'
********************************************

i=1
while(i<=32)

'set lon 100 130'
'set lat 20 50'
'set t 'i''
'set mpdset cnworld'
**'set mpdset cnriver '
**'set map 1 1 9'
'set xlopts  1 6 0.20'
'set ylopts  1 6 0.20'
'set mpdset cnworld'
'set parea 1 10 1 8'
'set xlopts  1 4 0.20'
'set ylopts  1 4 0.20'
'set clopts -1 -1 0.13'
'set xlab off'
'set ylab off'
'set grads off'
'set grid off'

'set gxout shaded'
'set cmin 10'
'set clevs      10 15 20 25 30 40 50 60 70 80 90 100 '
'set ccols 0     4 11 5  13 3  10 7  12 12  8  2  6 '
'cnbasemap  oacres(prc,r.2)'
'cbarn.gs'

'q time'
x=subwrd(result,3)
'draw title 'x' 6hr rain'
******坐标轴的设定***********************************************************
'run xylab3 100 130 5 20 50 5 0.15 1 5 0.10 0.10'
*****************
'print'
'c'
i=i+2
endwhile
'disable print'
'reinit'
;

注意：如果画图的范围包括这个中国区域，那输出gmf 文件没有问题；

如果输出的范围比全国的小，仍然用cnbasemap会在grads中图形显示正确，但是打开gmf格式时不能正常打开，提示“internal error:polyline is too long cannot continue drawing”：

这时需要把输出的类型改为png（'printim E:\rain\plot\'x'.gif white'），这样就不会有错误了。

我画出的图如下：

mofangbao

不错 真有心~

SKSY

哇、、、非常感谢啊。。这个实例很好的指导作用

曾记否

楼主辛苦了~~

七夕路

谢谢楼主 很实用

hoyizzz

怎么fortran运行出错啊

小龙虾

这个贴很好呢，谢谢楼主

怿燃

弄来试试，反正跟我要做的是一样的。

atomee

好东西哇  学习一下

daxuanfang

好东东，学习学习~