登录后查看更多精彩内容~

您需要 登录 才可以下载或查看，没有帐号？立即注册

x

本帖最后由 zly4814624 于 2016-10-27 10:39 编辑

小生不才，这段时间一直在处理FY2E的资料，因为接触NCL的时间不是很长，所以一路坎坷。

期间一直找寻各种方法，寻求各种帮助，也是弄的焦头烂额。

因为作为非nc、grid等类型的模式数据，这种二进制卫星数据在使用ncl处理起来本身就没有他们来的方便。

再加上卫星中心提供的这种数据本身没有经纬度的信息，而是单独形成了一个经纬度对照表，更是增加了一些难度。

好在。

功夫不负有心人。

我在众多好心人的帮助下。

初步完成了运用ncl处理FY2e的AWX格式数据的脚本。

情怀到此，话不多说。

运用到的数据：FY2E_SEC_IR1_LCN_20160917_1930.AWX，正式的数据的分辨率为512*512，但是AWX的卫星数据有一个特点，就是有头文件，这个数据的头文件的记录长度为6，所以头文件的字符占用个数为6*512，所以这个数据的总字符量为518*512。在读取数据的时候，必须要将头文件的字符给略去。除了这个数据，我们还需要再添加一个数据：land.dat，这是国家卫星中心提供的经纬度对照表，分辨率为512*512，正好与上面的数据完美配对。由于上面数据的特殊性，我们只能通过使用经纬度对照表来将经纬度信息加载进去。

begin

system("date")

configPath="$WORKDIR/config"

configName="awx.config"

inputfileName = getArgsPara(params,"FY2E_SEC_IR1_LCN_20160917_1930.AWX")

datafileName = getArgsPara(params,"land.dat")

;*************************************读取数据****************************

fl  =  getArgsPara(params,"inputfileName")

f   = fbindirread(fl,0,(/518,512/),"byte")               ；读取二进制数据文件

vaild = byte2flt(f(:,:))      ; [512] x [512]

tbb  = new((/512,512/),float)

do i=6,517

do j=0,511

tbb(i-6,j) = vaild(i,j)                                      ；跳过头文件数据，将正式数据存入数组

end do

end do

do i=0,511

do j=0,511

index=tbb(i,j)

if (index .lt. 0) then

tbb(i,j) = tbb(i,j)+255

else

tbb(i,j) = tbb(i,j)

end if

end do

end do

;********************************读取兰伯特经纬度对照表数据*************************

f2       = getArgsPara(params,"datafileName")

latlon2d = fbindirread(f2,-1,(/512,512,2/),"float")

dimlonlat=dimsizes(latlon2d)                              ;XLONG第一维时间，第二位纬度，第三维精度

nlat=dimlonlat(0)

mlon=dimlonlat(1)

lon2d  = latlon2d(:,:,0)

lat2d  = latlon2d(:,:,1)

;***********************************设置数据属性****************************

lat2d@units               = "degrees_north"

lon2d@units               = "degrees_east"

tbb@lat2d                 = lat2d

tbb@lon2d                 = lon2d

tbb@units                 = ""

tbb@long_name             = ""

tbb@coordinates           = "lat2d lon2d"

minvalue=min(tbb)

maxvalue=max(tbb)

tbb@_FillValue            = -999.

;************************************设置配置文件****************************

configFile= configPath+"/"+configName

;等经纬投影经纬度配置

truelat1 = stringtofloat(readParamInfo(configFile,"TRUELAT1"))

truelat2 = stringtofloat(readParamInfo(configFile,"TRUELAT2"))

stdlon   = stringtofloat(readParamInfo(configFile,"STDLON"))

; ; *************************************设置工作台****************************

res                                         = True

res@gsnDraw                                 = False

res@gsnFrame                                = False

res@gsnMaximize                             = True

res@gsnPaperOrientation                     = "portrait"

res@gsnAddCyclic                            = False

res@pmTickMarkDisplayMode                   = "Always"

;*************************************网格设置****************************

setMapProject(res,truelat1,truelat2,stdlon,lat2d,lon2d,mlon,nlat)

res@mpGridAndLimbOn                     = False

res@mpLimitMode                         = "LatLon"

res@mpMinLatF                           =   5.0

res@mpMaxLatF                           =  40.0

res@mpMinLonF                           =  89.01

res@mpMaxLonF                           = 129.025

res@trGridType                              = "TriangularMesh"

res@cnFillOn                                = True

res@cnFillMode                              = "RasterFill"

res@cnLinesOn                               = False

res@cnLineLabelsOn                          = False

; 坐标轴属性

res@tmXBOn                                     = False    ;坐标刻度隐藏

res@tmXTOn                                     = False

res@tmYLOn                                     = False

res@tmYROn                                     = False

res@tmXBBorderOn                            = False    ;四周边框

res@tmXTBorderOn                            = False

res@tmYLBorderOn                            = False

res@tmYRBorderOn                            = False

;色标显示

res@lbLabelBarOn                            = False

res@lbLabelAutoStride                       = False

; 垂直色标条

res@lbOrientation                           = "Vertical"

;自定义等值线

res@cnLevelSelectionMode                    = "ExplicitLevels"

; *************************************色标设置****************************

colors = (/(/255,255,255/),(/  0,  0,  0/),(/33,33,33/),(/68,65,68/),(/132,128,134/),(/201,192,201/),(/230,222,232/),\

(/235,166,195/),(/239,142,160/),(/244,118,123/),(/245,73,73/),(/220,79,69/),(/153,118,78/),(/181,145,80/),\

(/224,187,89/),(/239,202,92/),(/219,207,130/),(/205,212,165/),(/191,217,199/),(/186,218,209/),(/163,213,234/),\

(/156,203,228/),(/146,190,220/),(/144,187,217/),(/140,182,214/),(/132,172,208/),(/130,169,206/),(/116,149,193/),\

(/108,139,186/),(/103,132,182/),(/103,132,182/),(/95,121,175/),(/95,121,175/), (/84,106,165/),(/74,93,157/),\

(/68,86,152/),(/65,82,149/),(/59,73,143/),(/49,60,133/),(/45,55,129/), (/43,52,127/),(/36,45,121/),(/34,43,120/),\

(/35,42,118/),(/33,39,116/),(/27,32,110/),(/22,25,104/),(/19,21,101/),(/17,18,98/),(/13,13,94/),(/11,11,92/),\

(/10,11,91/),(/10,11,91/),(/10,11,91/),(/10,11,91/),(/10,10,91/),(/11,12,92/),(/12,12,93/),(/12,12,93/),\

(/11,11,92/),(/10,10,91/),(/10,10,91/)/)*1.0 ; we multiply by 1 to make colors float

res@cnLevels                                = fspan(50,250,59)

res@cnFillColors                            = ispan(61,2,1)

res@cnConstFLabelTextDirection              = "Down"

cmap = colors/255.             ; normalize (required by NCL)

;*************************************绘制设置****************************

;图片输出目录

outputPath="/home/suzhou/AWX"

checkOrCreateDir(outputPath)

;产品文件名

outputfileName="FY2E_SEC_IR1_LCN_20160917_1930"

;创建工作台

wks_type = "png"

wks_type@wkWidth                                = 1000          ;工作台宽度

wks_type@wkHeight                               = 1000        ;工作台高度

wks = gsn_open_wks(wks_type,outputPath+"/"+outputfileName)

contour = gsn_csm_contour_map(wks,tbb,res)

gsn_define_colormap(wks,cmap)

;叠加shp底图

mapres=True

addEastChinaShpMap(mapres)

chinamap = add_china_map(wks,contour,mapres)

draw(contour)

frame(wks)

delete([/wks,res,contour/])

system("date")

end

绘制出来的效果如下图，通过使用micaps直接载入AWX数据绘制的结果进行对比，可以达到100%相同的效果。

FY2E_SEC_IR1_LCN_201609171930_AWX.png (889.87 KB, 下载次数: 1)

2016-10-27 09:50 上传

总结：

1、该数据的难点在于二进制的数据类型的辨析，和头文件的占用个数，一般只要知道头文件一共占有多少字符量和正式数据的数据量，就可以成功的读取数据文件，运用geosat.exe可以辨析数据文件的分辨率，而对于头文件的信息，可以运用fortran程序进行读取，本论坛有大神提供过这种方法，水友们可以去找找。

2、由于AWX格式数据有多重分辨率和多种投影方式，所以国家卫星中心只提供了512*512这一种分辨率的经纬度对照表，那么对于其他分辨率的数据就不能再使用经纬度对照表了，需要自己手动创建经纬度了。

这个脚本比较简陋，有些地方可能存在需要改进的地方，水友们如果有不明白的地方或者有需要改进的地方，请多多提出意见。

另外，这种帖子可以申精么，第一次发这种帖子，嘿嘿。