请选择 进入手机版 | 继续访问电脑版
爱气象,爱气象家园! 

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 54595|回复: 62

[经验总结] 运用NCL处理风云卫星2E的AWX格式数据的总结

  [复制链接]

新浪微博达人勋

发表于 2016-10-27 10:12:10 | 显示全部楼层 |阅读模式

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

您需要 登录 才可以下载或查看,没有帐号?立即注册 新浪微博登陆

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

总结:
1、该数据的难点在于二进制的数据类型的辨析,和头文件的占用个数,一般只要知道头文件一共占有多少字符量和正式数据的数据量,就可以成功的读取数据文件,运用geosat.exe可以辨析数据文件的分辨率,而对于头文件的信息,可以运用fortran程序进行读取,本论坛有大神提供过这种方法,水友们可以去找找。
2、由于AWX格式数据有多重分辨率和多种投影方式,所以国家卫星中心只提供了512*512这一种分辨率的经纬度对照表,那么对于其他分辨率的数据就不能再使用经纬度对照表了,需要自己手动创建经纬度了。


这个脚本比较简陋,有些地方可能存在需要改进的地方,水友们如果有不明白的地方或者有需要改进的地方,请多多提出意见。
另外,这种帖子可以申精么,第一次发这种帖子,嘿嘿。

评分

参与人数 1金钱 +10 贡献 +2 收起 理由
river + 10 + 2 赞一个!

查看全部评分

本帖被以下淘专辑推荐:

密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-10-27 10:42:20 | 显示全部楼层
感谢分享,很棒的教程,连我这个不会NCL的都基本看懂了
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-10-27 10:45:59 | 显示全部楼层
river 发表于 2016-10-27 10:42
感谢分享,很棒的教程,连我这个不会NCL的都基本看懂了

哈哈,我也是从一个初学者的角度来解释脚本的。共同学习。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-10-27 11:33:29 | 显示全部楼层
加上省界就更好了
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-10-27 13:05:37 | 显示全部楼层
楼主的精神很伟大!感谢分享
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-10-27 13:05:54 | 显示全部楼层
正在学习NCL,这个对我启发很大~~~
密码修改失败请联系微信:mofangbao

新浪微博达人勋

 楼主| 发表于 2016-10-27 14:06:26 | 显示全部楼层

这个也有专门的属性语句的。
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-10-27 15:46:49 | 显示全部楼层
棒棒哒!!!支持!!!
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-10-27 15:51:09 | 显示全部楼层
感谢楼主的分享~
密码修改失败请联系微信:mofangbao

新浪微博达人勋

发表于 2016-10-28 11:24:47 | 显示全部楼层
6666666666666666666666666
密码修改失败请联系微信:mofangbao
您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

Copyright ©2011-2014 bbs.06climate.com All Rights Reserved.  Powered by Discuz! (京ICP-10201084)

本站信息均由会员发表,不代表气象家园立场,禁止在本站发表与国家法律相抵触言论

快速回复 返回顶部 返回列表