- 积分
- 11819
- 贡献
-
- 精华
- 在线时间
- 小时
- 注册时间
- 2013-4-26
- 最后登录
- 1970-1-1
|
登录后查看更多精彩内容~
您需要 登录 才可以下载或查看，没有帐号？立即注册
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%相同的效果。
总结：
1、该数据的难点在于二进制的数据类型的辨析，和头文件的占用个数，一般只要知道头文件一共占有多少字符量和正式数据的数据量，就可以成功的读取数据文件，运用geosat.exe可以辨析数据文件的分辨率，而对于头文件的信息，可以运用fortran程序进行读取，本论坛有大神提供过这种方法，水友们可以去找找。
2、由于AWX格式数据有多重分辨率和多种投影方式，所以国家卫星中心只提供了512*512这一种分辨率的经纬度对照表，那么对于其他分辨率的数据就不能再使用经纬度对照表了，需要自己手动创建经纬度了。
这个脚本比较简陋，有些地方可能存在需要改进的地方，水友们如果有不明白的地方或者有需要改进的地方，请多多提出意见。
另外，这种帖子可以申精么，第一次发这种帖子，嘿嘿。
|
评分
-
查看全部评分
|
[复制链接]
窥视卡
雷达卡
发表于 2016-10-27 10:12:10
提升卡
置顶卡
沉默卡
喧嚣卡
变色卡
抢沙发
千斤顶
显身卡
楼主
发表于 2016-10-27 11:33:29