日志
卫星tbb资料处理
热度 1 |
前段时间搜索了论坛中的tbb资料的处理方法，有很多的帖子，这里列出一个讲的比较完善的帖子http://bbs.06climate.com/forum.php?mod=viewthread&tid=4982。但是如果直接使用tbb变量画图，与实际是不符的，需要减去173.15k（得到摄氏温度）才与micaps查看的数值相当，但是搜索了所有的帖子也没有找到原因。
论坛中直接使用grads处理awx格式的tbb数据方法如下：
先根据卫星的型号设置正确的ctl（不同的卫星数据经纬度是不同的，这里以fy2e为例）：
dset J:\tbb\FY2E_TBB_IR1_OTG_20100530_0700.AWX
title 2010 twpice fy2e_tbb
options template
options yrev
* yrev表示y轴反向
undef 0.0
xdef 1201 linear 45 0.1
ydef 1201 linear -60 0.1
zdef 1 linear 1000 1
tdef 1 linear 23z30may2010 1hr
vars 1
tbb 1 -1,40,1 Geopotential Heights (m)
endvars
可以编写绘图gs，如下：
'open J:\tbb\tbb2010053007.ctl'
'set mpdset cnworld'
'set csmooth on'
'set grads off'
'set grid off'
'set lon 50 120'
'set lat 20 60'
'set gxout shaded'
****自定义阴影图颜色-------------------------------
'set rgb 21 249 252 215'
'set rgb 22 209 249 205'
'set rgb 24 167 250 200'
。。。。。有很多。。。。。
*'set rgb 255 10 10 10'
*'set rgb 256 0 0 0 '
****
'd tbb-173.15' ; ***参考了micaps数值显示后发现应该减去100而不是273.15，得到的是摄氏温度
'cbarn 0.8 1 10.6'
'printim J:\tbb\tbb2010053007h.emf x1024 y768 white'
后来请教了师姐，终于找到了少见100的原因，期间还发现了直接使用grads打开awx数据的问题，所以前面给的那个链接的帖子，虽然可以画出图来，但是其实与实际是有一个小小的偏移的，是不对的。
因为awx格式的卫星数据是由2个记录的数据头，而直接用grads
处理时却忽略了这两个记录，导致与实际情况有了一个小偏移。
(1)首先少见100的原因（减173.15而不是273.15得到摄氏温度）
由于awx中的是以单个字符来存储数据的，数据范围只能在0-255之间变化，实际开氏温度会超过这个温度，于是在存储的时候，
将原始的数值减去100然后再将其转化成对应的字符存储在awx中，这样用一个字节就可以代表原来需要2个字节的整型温度，
大幅度的减小了存储空间。之所以是100，是由于数据的一级数据头中的格点数据基准值设置为100。
(2)针对直接使用grads画tbb图像有漂移，是由于卫星数据文件头的原因导致的。所以是不能用grads代开awx格式数据然后直接
画图，必须首先对原始数据进行处理，将tbb数据值读取出来，然后再作图。
其读取程序如下：
character*1 ch(1201)
integer*4 tbb(1201,1201) !real tbb(1201,1201)
open(40,file='j:/tbb/tbb.grd',form='binary')
open(01,file=‘J:\tbb\FY2E_TBB_IR1_OTG_20100530_0700.AWX,form='binary',access='direct',recl=1201)
do i=1,1201
read(01,rec=i+2)ch
do j=1,1201
tbb(j,i)=ichar(ch(j))+100 !tbb(i,j)=ichar(ch(j))+100.0
enddo
enddo
write(40)(tbb(j,i),j=1,1201),i=1201,1,-1)
end
上述程序中，tbb即为黑体亮度温度，单位是K氏温度。
然后写ctl文件，这是ctl中不需要yrev ，在fortran中已经改变了存储方向了：
dset J:\tbb\FY2E_TBB_IR1_OTG_20100530_0700.AWX
title 2010 twpice fy2e_tbb
xdef 1201 linear 45 0.1
ydef 1201 linear -60 0.1
zdef 1 linear 1000 1
tdef 1 linear 15z30may2010 1hr
vars 1
tbb -1,40,4 99 tbb(K)
endvars
tbb 1 99 tbb(K)
还有如果是在fortran程序中定义的tbb是实数的话，那么相应的在ctl文件中用的是tbb 1 99 tbb(K)。具体的原因可以看看grads描述文件中对变量中参数的设置含义。如果是form='binary'的话，可以直接配合着recl=1201使用；但是如果form='unformatted'，并且编译器是ivf的话，由于它默认的是4个字节长度的倍数，即recl=1201，其实表示的是4804个字节；所以这个时候需要在fortran编译的时候改一下设置project--settings----fortran（category选fortran data)----data options中在Use Bytes as RECL= Unit for Unformatted Files 选项前打钩；这样设置以后recl=1201就表示是1201个字节，就不会fortran运行错误。当然如果简单起见，也可以直接配合 form='binary',recl=1201来运行。
还有一点事提醒注意的：卫星tbb资料的存储是从西往东从北往南，先存储纬向方向的再存经向。（第一条记录是60°N，经度分别为45°E，45.1°E，45.2°E，。。。，164.9°E，165°E；第二条记录是59.9°N，经度分别是45°E，45.1°E，45.2°E，。。。，164.9°E，165°E；。。。；最后一条记录是60°S，经度分别是45°E，45.1°E，45.2°E，。。。，164.9°E，165°E）。
以上是理解不知道有没有不妥之处，如有不对的地方希望指正。