利用micaps6小时降水资料和国家卫星中心的TBB画实况

clare

一般写技术总结或是论文的时候都需要画出实况降水，露珠在论坛前人的基础上，在美化方面下了一些功夫。特此把俺作图的所有步骤和内容
都分享给大家，以此表明我爱气象家园。
————————————————————————————————————————————————————————————————
因为本人叙述事情能力差，所以基本的micaps资料处理方法：请参见http://bbs.06climate.com/forum.php?mod=viewthread&tid=13648
————————————————————————————————————————————————————————————————
我这儿是将15-16日两天的降水求和后画图，数据是用excel计算得出（不会用fortran处理）
感兴趣的话，我这儿也给出我的处理步骤（大能可以忽略第二部分）：
————————————————————————————————————————————————
第一部分，用Fortran读数据生成grd,同时将所有数据读成rainning.txt，方便后面用excel处理：

program read_surface
          parameter(nt=36)
          integer::n,num,n1,n2,n3,n4,n5
           real,allocatable::lon(:),lat(:),rd(:),temp(:)    !动态数组分别用来储存站点的降水
           character*8,allocatable::sta(:)      !动态数组，储存站号，站号必须是8个字符，5个的话会stnmap出不出来   
        character*12 filename(nt)        !用于储存nt个时次的文件名
           open(1,file='filename.txt')
           do i=1,nt
               read(1,*) filename(i)
               print*,filename(i)
            enddo
            close(1)               
        open(2,file='rain6h.grd',form='binary')  !储存72个时次的降水
                open(9,file='rainning.txt')
       do k=1,nt
        open(3,file='d:\lunwen\rain\'//filename(k))
        do i=1,13
                read(3,*)
           end do
                read(3,*) n5,n !将该时次的站点数赋值于n
                print*,n
        allocate(lat(n))  !分配动态数组
        allocate(lon(n))
        allocate(rd(n))
        allocate(sta(n))
                allocate(temp(n))
               
        close(3)
        open(4,file='d:\lunwen\rain\'//filename(k))   !获知n后，重新读取数据。
        do i=1,14
                read(4,*)
            end do
          do i=1,n        !只读降水，别的要素暂略
             read(4,*) sta(i),lon(i),lat(i),temp(i),rd(i)
              tim=0.0;nlev=1;nflag=1
              write(2)  sta(i),lat(i),lon(i),tim,nlev,nflag,rd(i)  !一个时次的降水场输入完毕
                          write(9,200)  filename(k),sta(i),lat(i),lon(i),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)
                deallocate(temp)
        200 format(1x,A13,A10,1X,F5.2,1X,F8.2,1X,F8.2)
!释放动态数组
        close(4)
        enddo
         close(2)
                 close(9)
end program read_surface
——————————————————filename.txt

filename.txt

rainning.txt:

rainning.txt

———————————————————————————————————————————————————————————————
第二部分，用excel的数据透视表功能得出想要时段的总降水量。（露珠实在是说不清，所以具体的请百度一下，很简单的。附几张处理过程图）：
从rainning里面拷贝的想要的数据到excel里面：

某一日的降水

      
通过数据透视表处理后的：

excel处理后第二行即为总降水量

      
检验一下数据是否丢失：     

检查处理后结果，以防数据丢失

     
24rain.txt:
   

最终数据的格式

数据处理以后，用fortran转成grd格式：
program main
parameter(m=2401)
integer sta(m)
real jd(m),wd(m),a(m)
character aa
character*8 id(m)
!-----read data
open(1,file='d:\lunwen\rain\24rain.txt')
do i=1,m
        read(1,*)sta(i),wd(i),jd(i),a(i)
end do
close(1)

open (3,file='d:\lunwen\rain\zdh.txt')
open (4,file='d:\lunwen\rain\24rain.grd',form='binary')
!ccccccccccccccccccc
do j=1,m
id(j)=char(j)
tim=0.0
nlev=1
nflag=1
write(3,*)id(j),wd(j),jd(j),tim,nlev,nflag,a(j)
write(4)id(j),wd(j),jd(j),tim,nlev,nflag,a(j)
end do
tim=0.0
nlev=0
nflag=1
write(4)id(j-1),wd(j-1),jd(j-1),tim,nlev,nflag
close (3)
close (4)
end

——————————————————————————————————————————————————————————————————
第三部分，grid生成：
program gen_grid
!这里是定义的变量
integer,parameter::x=161,y=241,t=1
real grid
!程序开始
grid=0.25
open(1,file='rain24grid.grd',status='replace',form='binary')
do i=1,t
do j=1,y
        do k=1,x
                write(1) grid
        end do
end do
end do
close(1)
!程序结束
end

————————————————————————————————————————————————————————————————————
第四部分，画图（由于图中使用粉红斜线表示TBB小于32摄氏度的地方，必须使用grads2.1)。另外，由于前期没有把所有文件合并，所以造成gs中的多次调用。出图后懒得再弄，介意的小伙伴可以参照兰溪之水提供的方法

合并数据

，酱紫就不会太繁琐（注意时次别弄错就好）.
具体的gs:
'reinit'
*--------------------------------------
'open D:\lunwen\rain\1rain24h.ctl'
'open D:\lunwen\rain\1gridday.ctl'
*************基本底图设定**************************
'set mpdset cnhimap'
'set frame on'
'set mpdraw on'
'set poli on'
'set map 1 1 7'
'set csmooth on'
'set grid off'
'set grads off'
'set parea 1 10.5 1.0 8.0'
'set xlopts 1 7 0.15'
'set ylopts 1 7 0.15'
*************维数变量设定****************
'set lon 70 130'
'set lat 15 55'
'set t 1'
**************坐标设定***************************
'set xaxis 70 130 10'
'set yaxis 15 55 5'
*'set xlab off'
*'set ylab off'
*****************************************

'set gxout shaded'
'set cmin 0'
*'set clevs 1 3 5 7 9 11 13 15 17 19 21 23 25'
'set clevs  10 25 50 60 70 80 90 100 120 140 160 180 200 250'
'set rgb 20 255 255 255'
'set rgb 21 205 236 251'
'set rgb 22 150 212 243'
'set rgb 23 106 173 220'
'set rgb 24 72 148 181'
'set rgb 25 113 195 179'
'set rgb 26 73 168 113'
'set rgb 27 116 193 75'
'set rgb 28 206 219 86'
'set rgb 29 249 195 77'
'set rgb 30 246 130 51'
'set rgb 31 232 78 41'
'set rgb 32 215 39 40'
'set rgb 33 181 26 33'
'set rgb 34 146 21 25'

'set ccols 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 '

**'set xlpos -20'                                                                           
**'set ylpos -20'
'cnbasemap  oacres(g.2,r.1)'
'cbarm 0.8 0 2.5 0.1'
'close 2'
****************画TBB********************
'close 1'
'open D:\lunwen\awx\tbb.ctl'
'set mpdset cnworld cnriver '
*'set mproj lambert'  
'set frame on'
'set mpdraw on'
'set poli on'
'set map 1 1 7'
'set csmooth on'
'set grid off'
'set parea 1 10.5 1.0 8.0'
'set xlopts 1 7 0.15'
'set ylopts 1 7 0.15'
*------------------------------------
'set lon 70 130'
'set lat 15 55'
*------------------------------------
*'set xlint -10'
*'set ylint -10'
'set xlab off'
'set ylab off'

********必须grads2.1中用斜线或是点标示阴影********
i = 60
n = 60
cmd = 'set rbcols'
while (i<256)
  'set rgb 'n' 'i' 'i' 'i
  cmd = cmd%' '%n
  i = i + 15
  n = n + 1
endwhile
cmd

** define 2 colors
'set rgb 20 255 0 0'
'set rgb 21 215 170 198'

** define 2 patterns
'set tile 0 2 6 6 6 20'
'set tile 1 3 7 7 6 21'
'set rgb 22 tile 0'
'set rgb 23 tile 1'

'set gxout shaded'
'set cmax -62'
'set ccolor 5'
'set cthick 7'
'd ave(tbb,t=107,t=156)-173'
** overlay stippled red/blue shading of v greater than 5 and less than -5
'set clevs -32'
'set ccols 23 -1'
'd ave(tbb,t=107,t=156)-173'

*******************
*'lab_latlon'
*'cbarm 0.8 0 2.5 0.1'
'close 1'
*****************************************
'open D:\lunwen\will\fnll.ctl'
'set lon 70 130'  
'set lat 15 55'
'set lev 1000 10'
*'set t 1 9'
*如果有上面这句，会造成作图不在一张画面上
'set xlab off'
'set ylab off'
'set csmooth on'
*------------------------------------
'define h=HGTprs'
'define u=UGRDprs'
'define v=VGRDprs'
'define uu=ave(u,t=1,t=9)'
'define vv=ave(v,t=1,t=9)'
'define hh=ave(h,t=1,t=9)'
'define uv=ave(u,t=2,t=3)'
'define vu=ave(v,t=2,t=3)'
*------------------------------------
'set gxout contour'
'set clab forced'
'set lev 500'
'set cthick 7'
'set ccolor 11'
'set cint 4'
'd smth9(hh/10)'
*------------------------------------
'set gxout contour'
'set clab masked'
'set lev 500'
'set cthick 10'
'set ccolor 11'
'set clevs 588'
'd smth9(hh/10)'

'close 1'
*---------------画某时刻的风场---------------------
'open D:\lunwen\fnl1508.ctl'
'set lon 70 130'  
'set lat 15 55'
'set lev 850'
'set t 1'
'set xlab off'
'set ylab off'
'set csmooth on'

'define h=HGTprs'
'define u=UGRDprs'
'define v=VGRDprs'
'set gxout barb'
'set digsiz 0.04'
'set cthick 5'   
'set ccolor 1'   
'set lev 850'
'd skip(u*2.5,2);skip(v*2.5,2)'
*----------画南海-------
'basemap O 0'
'southsea_value -v'
'close 1'
*-------------------------------------------
*------神级坐标设置，如需使用请去掉*标注，前面也有哦------------------------------
*'run axis.gs -type b  -position o  -lfont 5 -lsize 0.1 -lthick 0.5 -interval 5 -sinterval 5 -suffix `3.'
*'run axis.gs -type L  -position o  -lfont 5 -lsize 0.1 -lthick 0.5 -interval 5 -sinterval 5 -suffix `3.'

'printim d:\lunwen\result\c5-26-32.png x1500 y1000 white'
;

_________________
说明：TBB的资料请参照http://bbs.06climate.com/forum.p ... 6&page=10#pid389105
gs中使用的cnhimap,请看http://bbs.06climate.com/forum.php?mod=viewthread&tid=28492
gs中使用的ctl如下：
                    

————————————————————————————————————————————————————————————————————
第五部分，其实我的作图有很多地方可以精简。譬如用excel处理数据那儿，相信fortran应该可以轻松做到；还有就是gs中因为没有提前cat所有micaps资料，造成画图的gs繁琐（因为无法用set t 来调控想要时次的数据）。俺读书少，所以真心不愿意再修改了。

传说中的谁

有图有真相，有过程，有源码，真是有理有据令人xing福。
感谢楼主分享，希望看到更多的原创帖子

river

感谢分享，楼主辛苦了。图很漂亮！

兰溪之水

挺好的，如果追求完满点，南海的小插图再修修~

有风吹过

楼主辛苦了，谢谢分享。

yanglei_nuist

正好要学习卫星，谢谢楼主了

去楼人空

傻蛋真厉害……

fujianbeige

楼主辛苦了，谢谢分享。

dcy_23

实用{:eb502:}{:eb502:}

shenqiankai

感谢楼主分享