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

气象家园

 找回密码
 立即注册

QQ登录

只需一步,快速开始

新浪微博登陆

只需一步, 快速开始

搜索
查看: 285|回复: 6

[源代码] cinrad 雷达数据3维显示

[复制链接]

新浪微博达人勋

发表于 7 天前 | 显示全部楼层 |阅读模式

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

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

x
本帖最后由 stream2011 于 2020-7-29 08:55 编辑

(1)针对cinrad 雷达数据常见的bin和bz2类型,进行检测判断并解析;
(2)对解析出来的反射率,速度值进行常规显示,如:
    show_lay(r)  # 分层显示反射率
    # show_3d(r, ei=1)  # 显示每层三维分布反射率
    # show_contour3d(r)  # 显示三维等值线反射率
    # show_mesh(r)  # 显示mesh 曲面
    # show_surf(r)  # 显示surf 曲面
    # show_point3(r)  # 显示所有层反射率值

主要模块:

class Radar:
    def __init__(self, path):
        self.path = path
        if path.endswith('.bz2'):
            newfilepath = path[0:len(path) - 4]
            with open(newfilepath, 'wb') as new_file, bz2.BZ2File(path, 'rb') as file:
                for data in iter(lambda: file.read(100 * 1024), b''):
                    new_file.write(data)
            file = open(newfilepath, 'rb')
        else:
            file = open(self.path, 'rb')

        file.seek(0)
        self.RawData = np.array([int(i) for i in file.read()])
        # self.Name = name[:-4]
        self.Count = int(len(self.RawData) / 2432)  # 数据体个数
        self.RawArray = self.RawData.reshape(self.Count, 2432)

        tta = np.array([[tvt[36],tvt[42],tvt[44], tvt[46], tvt[48],tvt[50],tvt[52],tvt[54],tvt[56],tvt[64],tvt[66],tvt[68],tvt[70],tvt[72],
                         tvt[37],tvt[43],tvt[45], tvt[47], tvt[49],tvt[51],tvt[53],tvt[55],tvt[57],tvt[65],tvt[67],tvt[69],tvt[71],tvt[73]]
                        for tvt in list(self.RawArray)])  # .astype(float)
        ttarr = tta[:, :14] + tta[:, 14:] * 256
        theta = ttarr[:, :2].astype(float)  # 方位角, 仰角
        theta = theta * (180/8/4096)
        ttarr = ttarr[:, 2:]
        # 方位角,仰角,层数,Reflect起始距离,Speed起始距离,Reflect库长,Speed库长,Reflect库数,Speed库数,Reflect数据位置指针,Speed数据位置指针,速度分辨率,vcp类型
        ele = ttarr[:, 0]  # .astype(int)  # 层数

        # 提取第一层数据
        self.refv = []  # 每层反射率值
        self.suv = []  # 每层速度值
        self.refxyz = []  # 每层反射率点坐标
        self.suxyz = []  # 每层风速点坐标
        self.eles = []  # 每层仰角
        for elei in np.unique(ele):
            li = np.array(np.where(ele == elei)).flatten()
            tvt1 = ttarr[li, :]  # 方位角—》(x,y)  仰角——>高度
            rawarray1 = self.RawArray[li, :]
            thetai = theta[li, :]  # 方位角, 仰角
            ref_xy = []
            ref_v = []
            su_xy = []
            su_v = []
            for thi in range(tvt1.shape[0]):  # 获取每条射线具体的值
                #  反射率
                if tvt1[thi, 5] > 0:
                    refd = (tvt1[thi, 1] + range(tvt1[thi, 5]) * tvt1[thi, 3])/1000  # 起始距离+ 库数*裤长 ,每个反射率值对应的距离,在根据方位角和仰角计算(x,y)
                    ref0 = rawarray1[thi, tvt1[thi, 7]: tvt1[thi, 7] + int(tvt1[thi, 5])]  # 反射率值
                    ref1 = (ref0.astype(float) - 2) / 2 - 32
                    ref1[ref1 < 0] = 0  # 计算结果< 0 的地方为0
                    ref1[((ref0 == 0) | ref0 == 1)] = 0  # 原位置为0或者1的为0
                    # ii = np.array(np.where(ref1 > 0)).flatten()
                    # ref1 = ref1[ii]
                    if len(ref_v) == 0:
                        ref_v = ref1
                    else:
                        ref_v = np.hstack((ref_v, ref1))

                    # ss = refd * np.sin(np.deg2rad(thetai[thi, 1])) + refd * refd / 17000
                    the = np.deg2rad(thetai[thi, 0])
                    if len(ref_xy) == 0:
                        costheta = np.cos(thetai[thi, 1])  # 仰角
                        refd2 = refd*refd/17000
                        ref_xy = np.vstack([refd * (np.cos(the)*costheta), refd * (np.sin(the)*costheta), refd * np.sin(np.deg2rad(thetai[thi, 1]))+refd2]).transpose()
                    else:
                        costheta = np.cos(thetai[thi, 1])  # 仰角
                        refd2 = refd * refd / 17000
                        ref_xy = np.vstack((ref_xy, np.vstack([refd *(np.cos(the)*costheta), refd * (np.sin(the)*costheta), refd * np.sin(np.deg2rad(thetai[thi, 1]))+refd2]).transpose()))  # 坐标

                #  速度
                if tvt1[thi, 6] > 0:  # 没有库数,则当前径向没有数据
                    sdd = (tvt1[thi, 2] + range(int(tvt1[thi, 6])) * tvt1[thi, 4])/1000  # 起始距离+ 库数*裤长 ,每个速度值对应的距离,在根据方位角和仰角计算(x,y)
                    sd0 = rawarray1[thi, tvt1[thi, 8]: tvt1[thi, 8] + int(tvt1[thi, 6])]  # 速度值
                    if int(tvt1[thi, 10]) == 2:  # 根据速度分辨率转化
                        sd1 = (sd0.astype(float) - 2) / 2 - 63.5
                        sd1[sd0 == 0] = 0
                    elif int(tvt1[thi, 10]) == 4:
                        sd1 = sd0 - 2 - 127
                        sd1[sd0 == 0] = 0
                    # ii = np.array(np.where(sd1 > 0)).flatten()
                    # sd1 = sd1[ii]
                    if len(su_v) == 0:
                        su_v = sd1
                    else:
                        su_v = np.hstack((su_v, sd1))

                    # ss = sdd * np.sin(np.deg2rad(thetai[thi, 1])) + sdd * sdd / 17000
                    if len(su_xy) == 0:
                        costheta = np.cos(thetai[thi, 1])  # 仰角
                        sdd2 = sdd * sdd / 17000
                        su_xy = np.vstack([sdd *(np.cos(the)*costheta), sdd * (np.sin(the)*costheta), sdd * np.sin(np.deg2rad(thetai[thi, 1]))+sdd2]).transpose()
                    else:
                        costheta = np.cos(thetai[thi, 1])  # 仰角
                        sdd2 = sdd * sdd / 17000
                        su_xy = np.vstack((su_xy, np.vstack([sdd * (np.cos(the)*costheta), sdd * (np.sin(the)*costheta), sdd * np.sin(np.deg2rad(thetai[thi, 1]))+sdd2]).transpose()))  # 坐标

                    #  谱宽
                    puk0 = rawarray1[thi, tvt1[thi, 9]: tvt1[thi, 9] + int(tvt1[thi, 6])]  # 谱宽值
                    puk1 = (puk0.astype(float) - 2) / 2 - 63.5
                    puk1[puk1 < 0] = 0  # 计算结果< 0 的地方为0
                    puk1[(puk0 == 0) | puk0 == 1] = 0  # 原位置为0或者1的为0
                    #  谱宽距离和速度一致
                    # pukd = sdd
            # 插值一层的数据
            if len(ref_xy) > 0:  # 显示反射率
                self.refv.append(ref_v)
                self.refxyz.append(ref_xy)
                self.eles.append(thetai[thi, 1])  # 每层仰角
                # tl = max(abs(ref_xy.min()), abs(ref_xy.max()))
                # grid_x, grid_y = np.mgrid[-tl:tl:200j, -tl:tl:200j]  # x方向在0-1上均匀生成200个数,y方向在0-1上均匀生成500个数
                # ref_z0 = griddata(ref_xy, ref_v, (grid_x, grid_y), method='nearest')
                # # grid_z1 = griddata(ref_xy, ref_v, (grid_x, grid_y), method='linear', fill_value=5)
                # # grid_z2 = griddata(ref_xy, ref_v, (grid_x, grid_y), method='cubic', fill_value=5)
                # plt.imshow(ref_z0, origin='lower')
                # plt.title('fanshelv')
                # plt.show()
            if len(su_xy) > 0:  # 显示速度
                self.suv.append(su_v)
                self.suxyz.append(su_xy)







`)6N]19T@ZO2H48IN0847LT.png

单层3维示意图

单层3维示意图

反射率第4层

反射率第4层
KML]XAB$@B5NWWG(1_H6GLA.png

反射率第一层

反射率第一层

所有层空间分布图

所有层空间分布图

反射率第3层

反射率第3层

反射率第二层

反射率第二层

反射率第九层

反射率第九层

cinrad_read_test.py

12.05 KB, 下载次数: 7, 下载积分: 金钱 -5

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

新浪微博达人勋

 楼主| 发表于 7 天前 | 显示全部楼层
待优化点:
(1)解析速度虽然比第一次优化了很多,但是赶c解析的还是慢很多;
(2)怎么将mlab三维展示结果能在web上进行交互展示,即可进行一些常规的旋转等操作;
(3)bz2文件现在是用生成文件进行跳转,虽然可以自动完成,但是如果通过文件流方式可能性能更好
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

发表于 7 天前 | 显示全部楼层
抢个沙发先,楼主真牛。谢谢楼主分享。怎么把三维的在web显示呢,
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 7 天前 | 显示全部楼层
xqls0203 发表于 2020-7-29 08:56
抢个沙发先,楼主真牛。谢谢楼主分享。怎么把三维的在web显示呢,

我也不晓得,主要是如果点多交互很容易出现卡顿
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 成长值: 0
0
早起挑战累计收入
发表于 前天 12:55 | 显示全部楼层
可以关注下 QuickEarth 气象数据三维可视化引擎
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 前天 16:13 | 显示全部楼层
mofangbao 发表于 2020-8-3 12:55
可以关注下 QuickEarth 气象数据三维可视化引擎

好的,谢谢,期待
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

新浪微博达人勋

 楼主| 发表于 前天 16:16 | 显示全部楼层
mofangbao 发表于 2020-8-3 12:55
可以关注下 QuickEarth 气象数据三维可视化引擎

另外请问下,通过雷达数据监测强对流区域,这个效果好不好???
密码修改失败请联系微信:mofangbao
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册 新浪微博登陆

本版积分规则

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

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

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