用NCL直接读取MICAPS4的资料并绘图

柿子柿子柿子

首先是丢一张图来展示结果。这个是某天某时某气压层的流场。刚出图的时候感觉还是不错的，但也被这密密麻麻的线恶心了一下，之后越看越觉得美感十足。后面用了稀疏（stMinDistanceF）发现不够这个有艺术感，就又去掉了。


接下来的话就是直接展示程序，懒得上传了，就直接分享出来了。也方便大家和自己以后的查询。

function readMicaps4Model(fileName:string)
begin

        f = cbinread(fileName, -1, "byte")

        discriminator  = tochar(f(0:3))            ; 合法数据关键字(mdfs)
        type              = toshort(f(4))              ; 数据类型(4/11)
        modelName   = tochar(f(6:25))           ; 模式名字
        element         = tochar(f(26:75))         ; 物理量
        description     = tochar(f(76:105))       ; 附加描述信息
        Extent           = tochar(f(178:277))      ; 扩展段(目前没有用)

        ; 层次
        cbinwrite("tmp.bin", f(106:109))
        level = cbinread("tmp.bin", -1, "float")
        system("rm tmp.bin")

        ; 年  月  日  起报时刻  时区  预报时效
        cbinwrite("tmp.bin", f(110:133))
        tmp = cbinread("tmp.bin", -1, "integer")
        system("rm tmp.bin")
        year          = tmp(0)        ; 年
        month       = tmp(1)        ; 月
        day           = tmp(2)        ; 日
        hour          = tmp(3)        ; 起报时刻
        timezone    = tmp(4)        ; 时区
        period        = tmp(5)        ; 预报时效
        delete(tmp)

        ; 经度
        cbinwrite("tmp.bin", f(134:145))
        tmp = cbinread("tmp.bin", -1, "float")
        system("rm tmp.bin")
        startLongitude        = tmp(0)  ; 起始经度
        endLongitude         = tmp(1)  ; 终止经度
        longitudeGridSpace = tmp(2)  ; 经度格距
        delete(tmp)
        ; 经度
        cbinwrite("tmp.bin", f(146:149))
        latitudeGridNumber = cbinread("tmp.bin", -1, "integer")        ; 纬向经线格点数
        system("rm tmp.bin")

        ; 纬度
        cbinwrite("tmp.bin",f(150:161))
        tmp = cbinread("tmp.bin", -1, "float")
        system("rm tmp.bin")
        startLatitude           = tmp(0)        ; 起始纬度
        endLatitude           = tmp(1)        ; 终止纬度
        latitudeGridSpace = tmp(2)        ; 纬度格距
        delete(tmp)
        ; 纬度
        cbinwrite("tmp.bin",f(162:165))
        longitudeGridNumber = cbinread("tmp.bin", -1, "integer")        ; 经向纬线格点数
        system("rm tmp.bin")

        ; 等值线
        cbinwrite("tmp.bin", f(166:177))
        tmp = cbinread("tmp.bin", -1, "float")
        system("rm tmp.bin")
        isolineStartValue = tmp(0)
        isolineEndValue  = tmp(1)
        isolineSpace       = tmp(2)
        delete(tmp)

        ; 坐标系
        Lat = fspan(startLatitude, endLatitude, longitudeGridNumber)
        Lon = fspan(startLongitude, endLongitude, latitudeGridNumber)
        Lat@long_name  = "Latitude"
        Lon@long_name = "Longitude"
        Lat@units   = "degrees-north"
        Lon@units  = "degrees-east"
        Lat!0          = "Lat"
        Lon!0         = "Lon"
        Lat&Lat      = Lat
        Lon&Lon     = Lon

        ; 数据区
        if (type .eq. 4) then
                ; 读数据
                cbinwrite("tmp.bin",f(278:))
                dat = cbinread("tmp.bin", (/longitudeGridNumber, latitudeGridNumber/), "float")
                system("rm tmp.bin")
                ; 设置信息
                dat!0     = "Lat"
                dat!1     = "Lon"
                dat&Lat  = Lat
                dat&Lon = Lon
        else if (type .eq. 11) then
                ; 读数据
                cbinwrite("tmp.bin",f(278:))
                dat = cbinread("tmp.bin", (/2, longitudeGridNumber, latitudeGridNumber/), "float")
                system("rm tmp.bin")
                ; 模和角度
                Norm = dat(0,:,:)
                Angle = dat(1,:,:)
                delete(dat)
                ; 计算uv
                u = new((/longitudeGridNumber, latitudeGridNumber/), float)
                v = new((/longitudeGridNumber, latitudeGridNumber/), float)
                do j = 0, longitudeGridNumber-1
                        do i = 0, latitudeGridNumber-1
                                rad   = Angle(j,i) * get_pi("float") / 180
                                u(j,i) = Norm(j,i) * cos(rad)
                                v(j,i) = Norm(j,i) * sin(rad)
                        end do
                end do
                delete([/Norm, Angle/])
                ; 填充数据
                dat            = new((/2, longitudeGridNumber, latitudeGridNumber/), float)
                dat(0,:,:) = u
                dat(1,:,:) = v
                ; 设置信息
                Direction            = (/"u", "v"/)
                dat!0                 = "Direction"
                dat!1                 = "Lat"
                dat!2                 = "Lon"
                dat&Direction     = Direction
                dat&Lat             = Lat
                dat&Lon            = Lon
                delete(Direction)
        end if
        end if

        ; 设置数据信息
        dat@type              = tostring(type)
        dat@modelName    = tostring(modelName)
        dat@element         = tostring(element)
        dat@description     = tostring(description)
        dat@level              = tostring(level)
        dat@year              = tostring(year)
        dat@month           = tostring(month)
        dat@day               = tostring(day)
        dat@hour              = tostring(hour)
        dat@timezone        = tostring(timezone)
        dat@period            = tostring(period)
        dat@isolineStartValue = tostring(isolineStartValue)
        dat@isolineEndValue  = tostring(isolineEndValue)
        dat@isolineSpace       = tostring(isolineSpace)
        
        ; 返回数组
        ; 标量(:,:)
        ; 矢量(0,:,:) => u
        ;       (1,:,:) => v
        return dat
end


上面这段是读取MICAPS4数据的。因为M4和M3数据差别很大，昨天无聊就把好久之前拖着的程序给弄了，也正好NCL可以立马就展示出来我提取的对不对。结果看来还是不错了。下面的这个的话就是绘图的部分了，我放在了另一个文件（mian.ncl）里面。

load "readMicaps4Model.ncl"

begin

        ;HGT = readMicaps4Model("dat/HGT.500.18050108.000")

        UV = readMicaps4Model("dat/UV.500.18050108.000")

        wks = gsn_open_wks("png", "test")
        ;res = True
        ;res@gsnDraw        = False
        ;res@gsnFrame      = False
        ;res@cnFillOn         = False
        ;res@cnLinesOn      = True
        ;res@gsnMaximize  = True
        ;res@gsnAddCyclic  = False

        ;res@mpMinLatF = min(HGT&Lat)
        ;res@mpMinLonF = min(HGT&Lon)
        ;res@mpMaxLatF = max(HGT&Lat)
        ;res@mpMaxLonF = max(HGT&Lon)

        ;plot = gsn_csm_contour_map(wks, HGT, res)

        stres = True
        stres@gsnDraw          = False
        stres@gsnFrame         = False
        stres@gsnMaximize     = True
        stres@gsnAddCyclic     = False
        ;stres@stMinDistanceF = 0.035
        stres@mpMinLatF = min(UV&Lat)
        stres@mpMinLonF = min(UV&Lon)
        stres@mpMaxLatF = max(UV&Lat)
        stres@mpMaxLonF = max(UV&Lon)
        ;stplot = gsn_csm_streamline(wks, UV(0,:,:), UV(1,:,:), stres)
        splot = gsn_csm_streamline_map_ce(wks, UV(0,:,:), UV(1,:,:), stres)

        ;overlay(plot,stplot)
        draw(splot)
        frame(wks)
end



最后说一下有个不知道什么原因出错的地方。
readMicaps4Model这个函数里面，在byte转换成short的那里，我查到的short是2个byte，然后往临时文件里面写出2byte再按照short读取的话会得到一个长度为2的short数组……有谁知道这个是本身的bug还是什么情况？至于我为什么要先写出到一个临时文件再读回来，是暂时没找到怎么从byte数组转换成integer、float这类数据的方法。

weinihou

平流层的萝卜

加点颜色渲染，就梵高再现了

雨恬魅影

太厉害啦

lc2016

逛论坛都能笑死我

weety9394

谢谢楼主，有人知道python怎么读这个数据吗？我直接从micaps4软件下载下来的

柿子柿子柿子

weety9394 发表于 2020-5-20 16:42
谢谢楼主，有人知道python怎么读这个数据吗？我直接从micaps4软件下载下来的

以byteArray的方式读入数据，再用struct.unpack()解码

hPa

想请教下，怎样用Python获取micaps4的数据啊？

柿子柿子柿子

hPa 发表于 2022-1-21 11:00
想请教下，怎样用Python获取micaps4的数据啊？

你在论坛里搜一下，我记得有这帖子