micaps多时次单变量批处理及其后grads画图

平流层的萝卜

最近想要用grads出micaps的能见度资料的图，时间是2013年1月份31天每天8次，地区是华北区域，资料是micaps地面plot文件。虽然对于高手来说思路步骤不复杂，但自己毕竟略菜，遇到了一些困难，所幸通过学习论坛的帖子基本都解决了。作为回报，将自己操作过程记录一下，分享交流。（ps：高手路过即可~~）
整个流程基本分为以下2个部分：
一、用fortran批量读取多时次micaps的资料文件（31*8=248个），写成一个包含多时次物理量场的grd。
二、配之以ctl，.map，格点场grd，最后gs出图。
虽然流程如此，但因为micaps文件有其特殊性，还是有几个需要注意的地方， 比如缺测值为9999，要选出华北的站点、剔除以外的站点资料等。下面一一讲到。
下图是要处理的248个micaps原始站点数据截图，
   
通过观察micaps数据格式，可以看到从上到下依次是 ：头文件说明，时次和该时次参与观测的站点数的说明，作为主干的若干站点号、经纬度和其物理量数据说明。作为能见度的数据在第18列，（具体哪些数据怎么排的，参见micaps文件格式）。
1、有了上面的认识，可以直接用fortran读取文件248个文件的能见度资料了：
这里我用了清风http://bbs.06climate.com/forum.php?mod=viewthread&tid=3079批量读取文件的先进方法。
首先，将所有文件名通过cmd命令输入到一个filename.txt文件里，再用fortran读取之。比之前自己通过观察micaps文件名规律再设计循环变量读取方便了太多！而且简洁易行，两步到位：一、转到要处理的目录，我的是E:\gdd\plot ，二、输入以下语句：dir /b *.000>filename.txt 。就会E:\gdd\plot下面生成一个filename.txt，里边按行储存了所有后缀名是.000的micaps文件名。如下图。

很方便！
下边的事情就好办了，fortran语句如下:（实现了华北站点挑选、剔除缺测数据）parameter(nt=248)
integer::n,num
real,allocatable::lon(:),lat(:),vi(:)        !动态数组分别用来储存站点的经纬度和能见度值
real:: temp(14)
character*8,allocatable::sta(:)                !动态数组，储存站号，站号必须是8个字符，5个的话会stnmap出不出来
character*8 ss(112) !112为华北112个站点，用于储存站点号         
character*12 filename(nt)        !用于储存248个时次的文件名
open(11,file='sample.txt')  ！这个sample.txt里边按行储存了华北站点站号。
read(11,*)
do i=1,112
        read(11,*) ss(i)
enddo
close(11)
open(11,file='filename.txt')
do i=1,nt
        read(11,*) filename(i)
        print*,filename(i)
enddo
close(11)               
open(12,file='vis.grd',form='binary')        !储存248个时次的能见度场
do k=1,nt
        print*,'文件数',k
        open(11,file=filename(k))
        read(11,*)
        read(11,*) n,n,n,n,n        ！将该时次的站点数赋值于n
        allocate(lat(n))
        allocate(lon(n))
        allocate(vi(n))
        allocate(sta(n))
        !print*,n
        close(11)
        open(11,file=filename(k))   ！获知n后，重新读取。
        read(11,*)
        read(11,*)
        do i=1,n        !只读能见度，风速暂略        此循环为该时次的站点循环
                read(11,*) sta(i),lon(i),lat(i),(temp(j),j=1,14),vi(i)
                if(vi(i)==9999) then !若为缺测站点，直接pass
                continue
                else        !不是缺测站点，接下来判断是否属于华北地区
                num=0
                do ii=1,112
                        if(sta(i)==ss(ii)) then
                        num=1
                        !print*,k,'   ',sta(i),'    ',ss(ii),'    ',num
                        endif
                enddo
                !print*,vi(i)
               
                if(num==1) then   ！如果该站点属于华北并未缺测，则写入grd中。
                        tim=0.0;nlev=1;nflag=1
                        write(12)  sta(i),lat(i),lon(i),tim,nlev,nflag,vi(i)
                endif
                endif
        enddo
nlev=0
write(12) sta(n),lat(n),lon(n),tim,nlev,nflag
        !一个时次的能见度场输入完毕
        deallocate(lat)
        deallocate(lon)
        deallocate(vi)
        deallocate(sta)
        close(11)
enddo
end
最后生成一个vis.grd。第一步完成。
开始第二步：
2、配置各种文件出图之。
     这里从简吧，值得注意的事情（当然论坛之前就有人提过）：      a 待插格点背景grd对应的ctl必须与站点资料（vis.grd）对应的ctl时间设置相同。在此都是tdef 248 linear 02z01jan2013 3hr
     b 为了方便在色标后添加物理单位，修改cbarn.gs文件。在气象家园整合版中，在cbarn.gs文件中的'draw string 'xp' 'yt' 'hi 和'draw string 'xr' 'yp' 'hi  后边分别添加km，即可。    贴上画图的gs：
'open E:\gdd\plot\vis.ctl'
'open E:\gdd\plot\grid.ctl'
'enable print E:\gdd\plot\micaps_vis.gmf'

i=1
while(i<=248)
'set grads off'
'set lon 100 130'
'set lat 20 50'
'set lev 850'

'set t 'i''
'set mpdset cnriver '
'set map 1 1 9'
'set xlopts 1 6 0.2'
'set ylopts 1 6 0.2'

'set gxout shaded'
'set clevs 0.1 0.3 1 10 15 20 25 30'
'cnbasemap  oacres(g.2,vi)'
'cbarn.gs'

'q time'
x=subwrd(result,3)
'draw title 'x''

'print'
'c'
i=i+1
endwhile
'disable print'
'reinit'
;
最后的出图结果比较粗糙（其实有了边界控制maskout之类的语句，也可以在gs里实现只显示华北地区的功能，只不过这些还在学习中,所以图会不定时改进）。如下图：

多的不贴了，可以看到这几个时次河北南部的污染还是相当严重的。
顺便贴上能见度的标准（百度）：
1.能见度20-30公里 能见度极好 视野清晰
2.能见度15-25公里 能见度好 视野较清晰
3.能见度10-20公里 能见度一般
4.能见度5-15公里 能见度较差 视野不清晰
5.能见度1-10公里 轻雾 能见度差 视野不清晰
6.能见度0.3-1公里 大雾 能见度很差
7.能见度小于0.3公里 重雾 能见度极差
8.能见度小于0.1公里 浓雾 能见度极差



PS：最后的gs和图在不定时改进中。







2014年4月11日更新：
在回答@F.Lee的问题的时候，顺便把自己以前的写grd的程序简化了一下，当时其实写的就不简洁，比较臃肿，后来也一直没把简化了的程序更新，对读者也会造成理解上的困难。
直接附上简化的基本思路，就不写代码了。


1,首先删除有关动态数组的声明及调用。因为sta(i),lat(i),lon(i),tim,nlev,flag,rd(i)完全没必要用数组，在获取站数的前提下，用站数循环，循环体内直接用变量sta,lat,lon,tim,nlev,flag,rd就好，省去了不少麻烦事。

2.无需打开两次filename(k) ，直接一次打开该文件，在里边先获取n站数，然后读取站点降水，写入grd即可。

mofangbao

终于出来了~好贴！讲的挺清楚

平流层的萝卜

mofangbao 发表于 2013-4-14 13:43

                               
登录/注册后可看大图

终于出来了~好贴！讲的挺清楚

谢清风~~！

kongfeng0824

好贴！！！！！！

曾记否

好好学习了一下，这个必须顶！楼主辛苦了！

river

很好很强大，学习了。我fortran学得不咋滴，想问个问题
allocate(lat(n))
        allocate(lon(n))
        allocate(vi(n))
        allocate(sta(n))
和 deallocate(lat)
        deallocate(lon)
        deallocate(vi)
        deallocate(sta)是什么意思啊，是用来运行剩下那些时次的？

平流层的萝卜

river 发表于 2013-4-14 19:10

                               
登录/注册后可看大图

很好很强大，学习了。我fortran学得不咋滴，想问个问题
allocate(lat(n))
        allocate(lon(n))

这个是fortran里边的动态数组，allocate(lat(n))的意思就是将lat数组的长度设置为n，deallocate（lat）就是将数组lat还原为未知，方便下一次读取站点资料继续用，因为每个时次站点数目不同，所以用动态数组比较好。

river

平流层的萝卜 发表于 2013-4-14 19:22

                               
登录/注册后可看大图

这个是fortran里边的动态数组，allocate(lat(n))的意思就是将lat数组的长度设置为n，deallocate（lat）就 ...

哦，原来是这么高级的东西啊。那就是write（12）就是已经把所有时次的资料都写进grd了，楼主为什么要注释   !一个时次的能见度场输入完毕。希望没有问的太水

平流层的萝卜

river 发表于 2013-4-14 19:32

                               
登录/注册后可看大图

哦，原来是这么高级的东西啊。那就是write（12）就是已经把所有时次的资料都写进grd了，楼主为什么要注释 ...

木有事，你可能没有看懂程序思路~~。外层循环是时次（do k=1,nt ），内层循环是站点（ do i=1,n ），因为micaps资料是按照时次而非站点储存。每到一个时次，打开一个micaps文件，同时输入该时次的若干个站点的能见度（称之为该时次的能见度场），因为不仅仅有一个时次，而是nt个，每循环一个时次，输入一个站点数n不一的能见度场，最后一共输入了nt个时次的能见度场。不知道我说的是否清楚？

river

平流层的萝卜 发表于 2013-4-14 20:20

                               
登录/注册后可看大图

木有事，你可能没有看懂程序思路~~。外层循环是时次（do k=1,nt ），内层循环是站点（ do i=1,n ），因为 ...

呵呵，明白了。只是你的这种实现方法没怎么用过，花费了一些时间。谢谢你的耐心解答