cinrad 雷达数据3维显示

stream2011 · 发表于 2020-7-29 08:44:22

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

您需要登录才可以下载或查看，没有帐号？立即注册

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)