NCL 绘制卫星AWX格式数据的问题

zly4814624

本人新手，运用ncl绘制风云卫星2E的AWX格式的数据，这是二进制数据，运用fbindirread读取，另外在国家卫星网上下载到了经纬度对照表，数据的像素大小是512*512的，分辨率是13公里/度。但是绘制出来的图非常不正常，如下图，看上去像是四张重复的，不知道哪里出了问题，希望有接触过的大神们帮我指点迷经，以下是我的代码：


load "$WORKDIR/include/library.ncl"


begin

        system("date")
        configPath="$WORKDIR/config"
        configName="awx.config"
        inputfileName = getArgsPara(params,"inputfileName")      ;读取AWX数据路径
        datafileName = getArgsPara(params,"datafileName")        ;读取经纬度对照表数据路径
        ;*************************************读取数据****************************
        fl  =  getArgsPara(params,"inputfileName")

        f   = fbindirread(fl,-1,(/512,512/),"integer")                     ;读取AWX数据
        tbb = int2flt(f(:,:))      ; [512] x [512]

        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)

        wks = gsn_open_wks("png","/home/suzhou/AWX")
        lat2d@units               = "degrees_north"
        lon2d@units               = "degrees_east"
        ;***********************************设置数据属性****************************
        tbb@lat2d                 = lat2d
        tbb@lon2d                 = lon2d
        tbb@units                 = ""
        tbb@long_name             = ""
        tbb@coordinates           = "lat2d lon2d"
        tbb@_FillValue            = -999.
        ; ************************************设置配置文件****************************
        configFile= configPath+"/"+configName
        ;等经纬投影经纬度配置
        truelat1 = stringtofloat(readParamInfo(configFile,"TRUELAT1"))
        truelat2 = stringtofloat(readParamInfo(configFile,"TRUELAT2"))
        stdlon   = stringtofloat(readParamInfo(configFile,"STDLON"))
        ;图片大小
        imgWidth=stringtofloat(readParamInfo(configFile,"IMGWIDTH"))
        imgHeight=stringtofloat(readParamInfo(configFile,"IMGHEIGHT"))
        ; ;*************************************设置工作台****************************
        wks_type = "png"
        wks_type@wkWidth                            = imgWidth          ;工作台宽度
        wks_type@wkHeight                           = imgHeight          ;工作台高度

        res = True

        res@gsnDraw                                 = False
        res@gsnFrame                                = False
        res@gsnMaximize                             = True
        res@gsnPaperOrientation                     = "portrait"
        res@gsnAddCyclic                            = False         ; regional data
        ;tfDoNDCOverlay should be set False (default) when lat2d and lon2d are used.
        ;res@tfDoNDCOverlay                         = True                        
        res@pmTickMarkDisplayMode                   = "Always"

        ;*************************************网格设置****************************
        setMapProject(res,truelat1,truelat2,stdlon,lat2d,lon2d,mlon,nlat)
        res@mpGridAndLimbOn                         = False

        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
        ;*************************************绘制设置****************************
        ;图片输出目录
        outputPath="/home/suzhou/AWX"
        checkOrCreateDir(outputPath)
        ;产品文件名
        outputfileName="FY2E_SEC_IR1_LCN_20160830_0430_AWX"
        ;创建工作台
        wks = gsn_open_wks(wks_type,outputPath+"/"+outputfileName)
        plot = gsn_csm_contour_map(wks,tbb,res)
        ; gsn_define_colormap(wks,cmap)
        ;叠加shp底图
        mapres=True
        addEastChinaShpMap(mapres)
        chinamap = add_china_map(wks,plot,mapres)

        draw(plot)
        frame(wks)

        myPrint("Save images:"+outputPath+"/"+outputfileName+"."+wks_type+" successfully!")
        system("date")
end

freekiller

我记得awx文件还有头信息，貌似没有看到读取哦

MeteoInfo

如果不是必须用ncl，可以用MeteoInfoLab，支持AWX格式，读取数据和绘图很简单。参考此贴：http://bbs.06climate.com/forum.p ... &extra=page%3D1

zly4814624

MeteoInfo 发表于 2016-9-24 08:33
如果不是必须用ncl，可以用MeteoInfoLab，支持AWX格式，读取数据和绘图很简单。参考此贴：http://bbs.06cli ...

谢谢您，但是必须要用NCL。

zly4814624

freekiller 发表于 2016-9-23 22:26
我记得awx文件还有头信息，貌似没有看到读取哦

是的是的。确实是有头文件，但是不知道怎么读取头文件，能指导一下么。

freekiller

zly4814624 发表于 2016-9-24 10:55
是的是的。确实是有头文件，但是不知道怎么读取头文件，能指导一下么。

可以找份文件说明文档，按照格式把头信息读取来。另外，论坛里有fortran，读取文件，把头信息扔掉，再输出文件。另外一种方法就是用ncl直接调用fortran, 读取画图一起处理。

zly4814624

freekiller 发表于 2016-9-24 18:03
可以找份文件说明文档，按照格式把头信息读取来。另外，论坛里有fortran，读取文件，把头信息扔掉，再输 ...

谢谢！。请问您说的文件说明文档是什么。

freekiller

zly4814624 发表于 2016-9-26 00:07
谢谢！。请问您说的文件说明文档是什么。

awx文件格式

zly4814624

freekiller 发表于 2016-9-26 09:20
awx文件格式

我设置了一段语句：
        aa   = fbindirread(fl,-1,(/6,512/),"byte")
        f     = fbindirread(fl,-1,(/512,512/),"byte")
        tbb  = byte2flt(f(:,:))      ; [512] x [512]

第一行的6是头文件的记录长度，第二行是读取数据，但是画出来的图像也不正常， 图像最顶端有一道多余的色条，那个好像是头文件信息被当做的数据也被绘制进去了，但是我已经设置了读取头文件了呀，请问不知道怎么回事。

zly4814624

zly4814624 发表于 2016-9-26 10:13
我设置了一段语句：
        aa   = fbindirread(fl,-1,(/6,512/),"byte")
        f     = fbindirrea ...

我通过新建数组略去了前6*512个字节的头文件，多余的色条问题解决了。
还有想请问一下，我想通过读取数据结构的方式来读取头文件信息，请问NCL能否实现这一行为。