使用NCL画风云卫星TBB数据

koala54

论坛里关于用grads画TBB（相当黑体亮度温度）数据的帖子已经有很多了，这里分享一些用NCL画TBB的经验，TBB数据来自国家卫星气象中心，下载方法可见如下帖子点我去看。为了方便验证做图结果，使用了clare的帖子新手分析卫星云图用grads画TBB图里面的例子。
作图思路：TBB是圆盘标称投影数据（何为圆盘标称投影？点我去看），投影坐标主要由中央经度和卫星高度两个参数确定，这样的数据分辨率一般是以实际距离m或km为单位，要用NCL画TBB主要问题就是把投影坐标转换到经纬度坐标。NCL可以直接打开hdf格式文件，用ncl_filedump查看一下文件，发现TBB数据在水平方向有2288*2288个格点，分辨率为5km（见文档1,P34）。这样的格点不是等经纬度的，可以通过格点参数计算出每个格点的经纬度（如用MeteoInfo的方法），或者根据卫星气象中心提供的经纬度查阅表(look-up table)读入格点的经纬度坐标。这里使用后一种方法，知道了格点的经纬度信息，就能直接用NCL画图了。

1. ;---------------------------------------------------
   
2. ; Plot FY2E TBB at 20130822_0000 (UTC)
   
3. ;---------------------------------------------------
   
4. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_code.ncl"
   
5. load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/gsn_csm.ncl"
   
6. 
   
7. begin
   
8. 
   
9. ;---- Read FY2 satellite data TBB.
   
10.         f         = addfile("FY2E_TBB_IR1_NOM_20130822_0000.hdf","r")
    
11.         tbb = f->FY2E_TBB_Hourly_Product
    
12. 
    
13.         LVR = tbb@LowerValidRange ; The valid range of data is [LVR,UVR].
    
14.         UVR = tbb@UpperValidRange
    
15. 
    
16. ;---- Read FY2 LatLon lookup table.
    
17.         f0                                 = "NOM_ITG_2288_2288(0E0N)_LE.dat"
    
18.         dims                         = (/2288,2288/)
    
19.         type                        = "float"
    
20.         ; Band 1 is lon: data1
    
21.         data1                        = fbindirread(f0,0,dims,type)
    
22.         ; Band 2 is lat: data2
    
23.         data2                         = fbindirread(f0,1,dims,type)
    
24.         
    
25.         lonCenter                 = 104.5                ; footpoint of FY-2E
    
26.         latCenter           = 0
    
27.         lon2d                         = data1+lonCenter
    
28.         lat2d                   = data2+latCenter
    
29.         lat2d@units         = "degrees_north"
    
30.         lon2d@units         = "degrees_east"
    
31. 
    
32. ;---- Add meta data for TBB.
    
33.         tbb@lat2d                 = lat2d
    
34.         tbb@lon2d                 = lon2d
    
35.         tbb@units                 = "degree Kelvin"
    
36.         tbb@long_name         = "Temperature of Bright Blackbody"
    
37.         tbb@coordinates = "lat2d lon2d"
    
38.         tbb@_FillValue  = -999.
    
39. 
    
40. ;---- Replace data outside the valid range with missing value.
    
41.         tbb_1d                                          = ndtooned(tbb)
    
42.         ind_notValid                         = ind(tbb_1d.lt.LVR.or.tbb_1d.gt.UVR)
    
43.         tbb_1d(ind_notValid)         = tbb@_FillValue
    
44.         tbb                                         = onedtond(tbb_1d,dimsizes(tbb))
    
45. 
    
46.         tbb = tbb-273.15        ;(convert degree Kelvin to degree Celsius)
    
47.         tbb@units = "degree Celsius"
    
48.         printMinMax(tbb,0)
    
49.         printVarSummary(tbb)
    
50. 
    
51. 
    
52. 
    
53. ;---- Begin to plot.
    
54.         wks = gsn_open_wks("png","TBB_plot")
    
55.         res = True
    
56. 
    
57.         res@gsnDraw                                 = False
    
58.         res@gsnFrame                                = False
    
59.         res@gsnMaximize                                = True
    
60.         res@gsnPaperOrientation                = "portrait"
    
61.         res@gsnAddCyclic                         = False         ; regional data
    
62.         ;tfDoNDCOverlay should be set False (default) when lat2d and lon2d are used.
    
63.         ;res@tfDoNDCOverlay                        = True                        
    
64.         res@pmTickMarkDisplayMode          = "Always"
    
65. 
    
66.         ;---- Map Set
    
67.         res@mpDataSetName                                = "Earth..4"
    
68.         res@mpDataBaseVersion                        = "MediumRes"
    
69.         res@mpOutlineSpecifiers                        = (/"China","China:Provinces"/)
    
70.         res@mpOutlineBoundarySets           = "NoBoundaries"
    
71.         res@mpNationalLineColor                        = "black"
    
72.         res@mpProvincialLineColor                = "black"
    
73.         res@mpGeophysicalLineColor          = "black"
    
74.         res@mpNationalLineThicknessF         = 3
    
75.         res@mpProvincialLineThicknessF  = 3
    
76.         res@mpGeophysicalLineThicknessF = 3
    
77. 
    
78.         res@mpMinLonF                                = 105
    
79.         res@mpMaxLonF                                = 128
    
80.         res@mpMinLatF                                = 18
    
81.         res@mpMaxLatF                                = 35
    
82. 
    
83.         res@trGridType                                = "TriangularMesh"
    
84.         res@cnFillOn                                 = True
    
85.         res@cnFillMode                                 = "RasterFill"
    
86.         res@cnLinesOn                                 = False
    
87. 
    
88.         res@gsnLeftString                        = "TBB"
    
89.         res@gsnRightString                        = "~S~o~N~C"
    
90. 
    
91.         ;---- Color set
    
92.         cmap  = read_colormap_file("precip3_16lev")
    
93.         cmap1 = cmap(4::,:)
    
94.         cmap1(0,:) = cmap(0,:)
    
95.         cmap1 = cmap1(::-1,:)
    
96.         res@cnFillPalette = cmap1
    
97.         
    
98.         res@cnLevelSelectionMode = "ExplicitLevels"
    
99.         res@cnLevels = (/-30,-40,-50,-60,-70,-80/)
    
100.         res@lbOrientation = "vertical"
     
101. 
     
102.         plot = gsn_csm_contour_map_ce(wks,tbb,res)
     
103.         
     
104.         draw(plot)
     
105.         frame(wks)
     
106. 
     
107. end
     

复制代码

2016-05-15：
几点心得：
1）经纬度查阅表提供的是二维经纬度信息 ，即经度和纬度都有2288*2288个格点，在NCL中二维经纬度不能作为变量的坐标变量，而要放在变量的属性中（如 tbb@lat2d=lat2d，而不是tbb!lat2d=lat2d）;
2）知道二维经纬度后可以把图画在任意的地图投影里（脚本以极射赤面投影为例）；
3）经纬度查阅表中的格点是默认以0°E和0°N为起点，使用时要加上卫星观测的中央经纬度（图1）；
4）tbb变量的属性中有两个属性LowerValidRange和UpperValidRange分别代表有效数据的范围，范围之外的数据在脚本中设为了缺省值；
5）我看很多用grads画tbb的帖子中，最后把tbb-173和tbb+100分别转换成℃和K为单位；而我画hdf格式的tbb时，觉得tbb单位是K, tbb-273.15单位为℃时画出来结果和用grads画出结果比较一致，我没有找到相关的说明文档，大家知道的话请告诉我一下

2016-05-27 更新：
1）使用了兰溪之水版主重制的中国地图，边界准确，绘图快速，使用方便，在这里推荐给大家（去看看）；
2）地图投影更改为低纬地区常用的等距圆柱投影；
3）关于GrADS绘制风云卫星AWX格式的数据，为什么要将TBB-173和TBB+100分别转换成℃和K为单位，我的理解是：
    AWS文件是二进制格式，数据格点的空间分辨率为0.1°x0.1°，东西和南北方向各有1201个格点，空间范围为：45-165°E，-60°S-60°N。数据开头有两行说明，所以数据可以看作有1203行和1201列。为了节约存储空间，TBB数据并没有以整型（4字节）存储，而是以ASCII字符（1字节）存储，并以ASCII字符的标号表示TBB。因为ASCII字符标号的范围只有0-255，而TBB（单位K）的范围要超过这个范围，所以TBB先减去了100再将其转换成ASCII字符存储。如果直接用GrADS画AWX格式数据，数据会向北偏移0.2°（有两行说明的缘故），可以先用Fortran提取出TBB数据再画图（程序见附件）。
    HDF格式的TBB直接是用K为单位存储的，故可直接使用。
    提取AWS二进制格式卫星数据的Fortran程序：
   

1. program AWX2GRD
   
2. 
   
3.     implicit none   
   
4.     integer,parameter :: N=1201
   
5. 
   
6.     character*1 :: ch(N,N+2)
   
7.     integer :: tbb(N,N),i,j
   
8. 
   
9.     open(10,file='F:\TBB\FY2E_TBB_IR1_OTG_20120816_0000.AWX',form='unformatted',access='direct',recl=1201*1023)
   
10.     open(20,file='F:\TBB\FY2E_TBB_IR1_OTG_20100818_0000.grd',form='unformatted')
    
11.       
    
12.     read(10,rec=1) ch
    
13. 
    
14.     tbb = ichar(ch(:,3:N+2))
    
15.     tbb = tbb+100
    
16. 
    
17.     write(20)((tbb(j,i),j=1,N),i=N,1,-1)
    
18. 
    
19. end program AWX2GRD

复制代码

sgfzl

FY-2E自2015年7月1日起漂移至东经86.5度赤道上空，所以2015年7月1日以后的数据中心经度都要改成东经86.5度

yujie

谢楼主分享

曦尘

很好，谢谢分享！

锦叶seven

请问楼主，这里加载的load "$GEODIAG_ROOT/geodiag.ncl"是什么呀

freekiller

谢谢分享。      

koala54

锦叶seven 发表于 2016-5-16 14:04
请问楼主，这里加载的load "$GEODIAG_ROOT/geodiag.ncl"是什么呀

加载的是大气所董理老师制作的NCL地理数据诊断包，这里用来绘制地图底图（;----Map set后的内容），具体介绍可见http://www.lasg.ac.cn/Sylm/2012/2/ni52mtkrni.htm，下载及安装见https://github.com/dongli/geodiag

锦叶seven

koala54 发表于 2016-5-16 18:25
加载的是大气所董理老师制作的NCL地理数据诊断包，这里用来绘制地图底图（;----Map set后的内容），具体 ...

好的，谢谢楼主

狼牙李

感谢楼主分享

owen2121

真给力的教材啊  点赞

tiantianyixin

koala54 发表于 2016-5-16 18:25
加载的是大气所董理老师制作的NCL地理数据诊断包，这里用来绘制地图底图（;----Map set后的内容），具体 ...

你好，我想下载$GEODIAG_ROOT/geodiag.ncl脚本，但在官网上下载不成功，能否分享一下，谢谢你