X波段双极化相控阵天气雷达基数据的python读取方法
雷达数据的解析是雷达数据分析应用的基础，前段时间《新一代天气雷达》的作者张深寿老师推送了一篇天气雷达基数据解析的文章《Python读取和显示天气雷达基数据标准格式雷达基数据》，大家都受益匪浅，随着国内双极化相控阵雷达的应用兴起，有许多小伙伴希望我们能够科普双极化相控阵雷达的基数据解析，经过一段时间的整理，现在提供给各位小伙伴。
一、常用天气雷达基数据处理工具
随着气象雷达的发展，气象雷达数据处理领域产生了若干优秀的开源工具，雷达基数据格式也是百花齐放。
pyART 项目
Py-ART是一款Python工具包，里面包含了一组天气雷达算法和应用程序，Py-ART开始是大气遥感气候研究所（ARM）用来处理降雨和云雷达的数据的工具，现在已经被广泛应用于多种雷达和大气领域的数据处理、分析和检验。
其IO模块支持的格式：MDV二进制，NEXRAD-II二进制，Cfradial NETCDF格式，Sigmet二进制格式。
WRadLib 项目
WRadLib是一款开源的天气雷达数据处理程序，由Python语言编写。WRadLib可以帮助用户快速解决很多雷达数据处理工作，例如：数据读取、读取通用数据格式、坐标转换、降雨估计、杂波处理、衰减订正、数据可视化。
其IO模块支持的格式：Sigmet-IRIS二进制格式，OPERA-HDF格式，Radola(Radar-Online-Aneichung）德国气象雷达拼图格式。
pyCWR 项目
pyCWR是南信大大气物理学院专业雷达团队基于pyART开发的优秀国产气象雷达开源项目。
其IO模块支持的格式：metstar双偏振SAB，SA/SB/CB，SC/CD，CC/CCJ。
二、常用天气雷达的数据格式的结构
传统天气雷达的基数据一般都包括头信息和数据块，有些格式有1个总头信息块（例如：WSR88D格式），有些格式每个径向都有独立的径向信息快（例如：SA/SB和CB格式）。纳睿雷达（naruida）双极化相控阵天气雷达基数据格式采用了没有总头信息块，每个径向信息块包含全部信息的方式，目的是为了增强网络数据传输过程中冗余性。
头信息和数据快组成二进制文件，压缩方式包括：全局压缩和内部数据块压缩两种方式，全局压缩的优缺点都是必须全部解压后才能获得完整的头信息。内部数据快压缩方式的基数据文件常采用HDF和NETCDF格式的成熟文件存储协议。
MDV二进制全局头格式（美国）
MDV（Meteorological Data Volume）格式拥有全局头信息，特点是：可在文件末尾增加扩展数据块（chunk data），例如，回波的CAPPI-2d数据，可以在chuck上存放每个网格点对应的插值仰角，距站点的距离等，这些信息可用来计算雷达拼图或者其他产品。
Cfradial/CF2-NETCDF格式（美国）
针对极坐标采样探测系统设计，对微波和激光探测数据采用统一的数据格式，针对数值预报和同化系统做适应性开发，VDRAS同化系统采用这种雷达数据格式作为输入。
Sigmet-IRIS二进制全局头格式（欧美）
芬兰维萨拉公司气象雷达基数据格式，拥有独立信息块，每个径向拥有独立冗余的头信息，南京大学C-pole和珠澳S波段雷达采用这个格式。
NEXRAD-II / WSR98D /CINRAD-双极化SAB二进制全局头格式
美国NEXRAD组网S波段双极化雷达，中国CINRAD组网升级双极化SAB型，拥有全局头信息块，每个径向拥有独立的径向头信息。
CINRAD-SC/CD/CC二进制全局头格式
锦江SC/CD型，四创CC型，拥有全局头信息块，每个径向拥有独立的径向头信息。
CINRAD-SA/SB/CB/SC2.0和naruida-NRD二进制径向头格式
Metstar SA/B和CB型，无全局头信息块，每个径向有自己的径向头信息块。
三、naruida双极化相控阵天气雷达基数据格式的特点
纳睿雷达双极化相控阵天气雷达的基数据协议参考了欧美雷达基数据格式和我国CINRAD雷达网的6 种基数据格式，考虑到综合一体化气象观测的业务需求，结合当下分布式存储系统和数据并发传输的技术现状，和雷达系统未来网络化，多功能，自适应，软件化，智能化的发展趋势，基于自身VRHI扫描模式配置丰富多样的特点，采用了和国家气象局大气探测中心雷达SAB多普勒雷达基数据一致的数据组织方式，是一种适用于海量数据实时网络传输基数据协议，优点是：无全局数据头（去中心化），以径向探测Radial为最小检索单元，逐RHI扫描添加。径向2000库和8个变量的数组，无压缩100MB，体积可以比照SA/SB多普勒格式的3要素二进制文件40MB，SA/SB双偏振格式10要素100MB。
四、naruida双极化相控阵天气雷达基数据python读取
（1）python版本
3.6
（2）安装
git clone https://github.com/YvZheng/pycwr.git
sudo apt-get install python-numpy
pip install Cython
（3）数据读取代码
脚本位置：pycwr_master/pycwr/io/AXPT_NRD.py
from .BaseDataProtocol.AXPT_NRD_Protocol import dtype_AXPT0164
fillvalue_undef = -988.0
fillvalue_no_obs = -999.0
class AXPT0164BaseData(object):
"""
解码 Naruida双极化相控阵雷达的数据格式
"""
def __init__(self, filename, station_lon=None, station_lat=None, station_alt=None, debug=False):
"""
:param filename: radar basedata filename
:param station_lon: radar station longitude //units: degree east
:param station_lat: radar station latitude //units: degree north
:param station_alt: radar station altitude //units: meters
"""
#time0 = time.time()
super(AXPT0164BaseData, self).__init__()
self.debug = debug
self.filename = filename
self.station_lon = station_lon
self.station_lat = station_lat
self.station_alt = station_alt
self.fid = _prepare_for_read(self.filename) ##对压缩的文件进行解码
self._check_standard_basedata() ##确定文件是standard文件
#print(time.time() - time0) #0.05s
#time0 = time.time()
fixed_buf = self.fid.read(dtype_AXPT0164.RadarInfoHeaderBlockSize + dtype_AXPT0164.TaskConfigurationBlockSize) ##读取第一个固定头的信息
self.fid.seek(0, 0)
self.header = self._parse_BaseDataHeader(fixed_buf)
#print(time.time() - time0) ##0.0001秒
#time0 = time.time()
self.radial = self._parse_radials() ##维度次序==[仰角，方位]
#print(f' _parse_radials: {time.time()-time0} seconds') ##mask从c计算7秒
#time0 = time.time()
self.nrays = len(self.radial)
status = np.array([istatus['RadialState'] for istatus in self.radial[:]])
self.sweep_start_ray_index = np.where((status == 0) | (status == 3))[0]
self.sweep_end_ray_index = np.where((status == 2) | (status == 4))[0]
self.nsweeps = len(self.sweep_start_ray_index)
self.fid.close()
#print(time.time() - time0) ##0.003秒
def _check_standard_basedata(self):
"""
:param fid: file fid
:return:
"""
assert struct.Struct('H').unpack_from(self.fid.read(2), 0)[0] == 15, 'file in not a stardand AXPT0164 file!'
self.fid.seek(0, 0)
return
def _parse_BaseDataHeader(self, fixed_buf):
BaseDataHeader = {}
BaseDataHeader['RadarInfoHeaderBlock'], _ = _unpack_from_buf(fixed_buf, \
dtype_AXPT0164.RadarInfoHeaderBlockPos,
dtype_AXPT0164.BaseDataHeader['RadarInfoHeaderBlock'])
BaseDataHeader['TaskConfig'], _ = _unpack_from_buf(fixed_buf,
dtype_AXPT0164.TaskConfigurationBlockPos,
dtype_AXPT0164.BaseDataHeader['TaskConfigurationBlock'])
BaseDataHeader['RadarInfoHeaderBlock']['Longitude'] = BaseDataHeader['RadarInfoHeaderBlock']['Longitude'] * (360.0 / 65535.0)
BaseDataHeader['RadarInfoHeaderBlock']['Latitude'] = BaseDataHeader['RadarInfoHeaderBlock']['Latitude'] * (360.0 / 65535.0)
BaseDataHeader['RadarInfoHeaderBlock']['BeamWidthVert'] = BaseDataHeader['RadarInfoHeaderBlock']['BeamWidthVert'] * (10.0 / 65535.0)
BaseDataHeader['RadarInfoHeaderBlock']['BeamWidthHori'] = BaseDataHeader['RadarInfoHeaderBlock']['BeamWidthHori'] * (10.0 / 65535.0)
BaseDataHeader['RadarInfoHeaderBlock']['Height'] = BaseDataHeader['RadarInfoHeaderBlock']['Height'] * (1000.0 / 65535.0)
BaseDataHeader['RadarInfoHeaderBlock']['SiteName'] = BaseDataHeader['RadarInfoHeaderBlock']['SiteName']
BaseDataHeader['RadarInfoHeaderBlock']['Frequency'] = 9500 ##9500 MHz = 10e6 Hz
BaseDataHeader['TaskConfig']['ScanUtcTime'] = datetime.datetime.fromtimestamp(BaseDataHeader['TaskConfig']['VolumeStartTime'])
BaseDataHeader['TaskConfig']['Azimuth'] = BaseDataHeader['TaskConfig']['Azimuth'] * (360.0 / 65535.0)
BaseDataHeader['TaskConfig']['Elevation'] = BaseDataHeader['TaskConfig']['Elevation'] * (360.0 / 65535.0)
BaseDataHeader['TaskConfig']['StartGateRange'] = BaseDataHeader['TaskConfig']['StartGateRange'] * (1000.0 / 65535.0)
BaseDataHeader['TaskConfig']['UnitGateLength'] = BaseDataHeader['TaskConfig']['UnitGateLength'] * (1000.0 / 65535.0)
BaseDataHeader['TaskConfig']['ScanType'] = 0 ##多层RHI=5， 体扫=0. 实际扫描方式为v-RHI，写入时伪装为S型雷达的扫描方式v—PPI
BaseDataHeader['TaskConfig']['PulseWidth'] = 20 ##nanosec=10e-9sec
BaseDataHeader['TaskConfig']['NyquistSpeed'] = 26.0 ##26m/s
BaseDataHeader['TaskConfig']['NyquistDistance'] = BaseDataHeader['TaskConfig']['NyquistDistance'] * 1000.0 ##meter
self.AzimuthNumber = BaseDataHeader['TaskConfig']['AzimuthNumber']
self.Azimuth_list = []
self.ElevationNumber = BaseDataHeader['TaskConfig']['ElevationNumber']
self.Elevation_list = []
self.LengthOfGate = BaseDataHeader['TaskConfig']['LengthOfGate']
return BaseDataHeader
def _parse_radials(self):
radial = []
#time0 = time.time()
size_bytes = dtype_AXPT0164.RadialDataSize * self.AzimuthNumber * self.ElevationNumber ##无法少读取 RHI-v
base_data = self.fid.read(size_bytes)
#print(type(base_data), len(base_data)) ##77299200 = 16104* self.AzimuthNumber * self.ElevationNumber
#print(time.time() - time0)
Azimuth_list = []
Elevation_list = []
NyquistDistance_list = []
for el in np.arange(0, self.ElevationNumber):
#print(el)
for az in np.arange(0, self.AzimuthNumber):
RadialDict = {}
###
factor = az*self.ElevationNumber*dtype_AXPT0164.RadialDataSize + el*dtype_AXPT0164.RadialDataSize
buf_radial = base_data[factor : factor+dtype_AXPT0164.RadialDataSize]
###
buf_header = buf_radial[0 : dtype_AXPT0164.RadarInfoHeaderBlockSize + dtype_AXPT0164.TaskConfigurationBlockSize] ##固定头的信息
radial_header = self._parse_BaseDataHeader(buf_header)
#print(round(radial_header['TaskConfig']['Azimuth'],2), end=", ")
#print(round(radial_header['TaskConfig']['Elevation'],2), end=", ")
#print(radial_header['TaskConfig']['StartGateRange'])
RadialDict['Azimuth'] = radial_header['TaskConfig']['Azimuth']
if el == 0:
Azimuth_list.append(RadialDict['Azimuth'])
RadialDict['Elevation'] = radial_header['TaskConfig']['Elevation']
if az == 0:
Elevation_list.append(RadialDict['Elevation'])
NyquistDistance_list.append(radial_header['TaskConfig']['NyquistDistance'])
RadialDict['Seconds'] = radial_header['TaskConfig']['VolumeStartTime']
RadialDict['MicroSeconds'] = 0.0
###
if az == 0:
if el == 0:
RadialDict['RadialState'] = 3 #体扫开始
else:
RadialDict['RadialState'] = 0 #某仰角开始
elif az == self.AzimuthNumber - 1:
if el == self.ElevationNumber - 1:
RadialDict['RadialState'] = 4 #体扫结束
else:
RadialDict['RadialState'] = 2 #某仰角结束
else:
RadialDict['RadialState'] = 1
if self.debug:
if RadialDict['RadialState'] in [0, 3]:
ptime = pendulum.from_timestamp(RadialDict['Seconds'])
print(' ele=', int(RadialDict['Elevation']), ptime, ' az=', RadialDict['Azimuth'])
###
buf_data = buf_radial[dtype_AXPT0164.DataBlockPos : dtype_AXPT0164.RadialDataSize - 4]
RadialDict['fields'], RadialDict['mask'] = self._parse_single_radial_field_v1(buf_data)
radial.append(RadialDict)
continue
continue
self.Azimuth_list = np.array(Azimuth_list)
self.Elevation_list = np.array(Elevation_list)
self.NyquistSpeed_list = np.array([self.header['TaskConfig']['NyquistSpeed']] * self.ElevationNumber)
self.NyquistDistance_list = np.array(NyquistDistance_list)
return radial
def _parse_single_radial_field_v1(self, buf_data):
radial_var = struct.unpack(dtype_AXPT0164.MomentNum*dtype_AXPT0164.n_variable*'B', buf_data) ##0.0005s
radial_var = np.array(radial_var).reshape(dtype_AXPT0164.n_variable, dtype_AXPT0164.MomentNum)
radial_var = dict(zip([i[0] for i in dtype_AXPT0164.BaseDataHeader['RadialFieldDataBlock']], radial_var)) ##0.005s
radial_var['dBT'] = radial_var['dBT'][:self.LengthOfGate] * 106. / 255. -20. ##dBZ 86
radial_var['dBZ'] = radial_var['dBZ'][:self.LengthOfGate] * 106. / 255. -20. ##dBZ 86
radial_var['V'] = radial_var['V'][:self.LengthOfGate] * 101. / 255. -50. ##m/s 51
radial_var['W'] = radial_var['W'][:self.LengthOfGate] * 101. / 255. -50. ##m/s 51
radial_var['ZDR'] = radial_var['ZDR'][:self.LengthOfGate] * 16. / 255. -5. ##dB 11
radial_var['PhiDP'] = radial_var['PhiDP'][:self.LengthOfGate] * 181. / 255. ##deg 108 !
radial_var['KDP'] = radial_var['KDP'][:self.LengthOfGate] * 31. / 255. -10. ##deg/km 21
radial_var['CC'] = radial_var['CC'][:self.LengthOfGate] * 1.1 / 255. ##1 1.1 !
##
mask = np.array([False] * self.LengthOfGate)
mask[radial_var['dBT']==86] = True
##
radial_var['dBT'][radial_var['dBT']==86] = fillvalue_undef
radial_var['dBZ'][radial_var['dBZ']==86] = fillvalue_undef
radial_var['V'][radial_var['V']==51] = fillvalue_undef
radial_var['W'][radial_var['W']==51] = fillvalue_undef
radial_var['ZDR'][radial_var['ZDR']==11] = fillvalue_undef
radial_var['PhiDP'][radial_var['PhiDP']==181] = fillvalue_undef
radial_var['KDP'][radial_var['KDP']==21] = fillvalue_undef
radial_var['CC'][radial_var['CC']==1.1] = fillvalue_undef ##0.001s
return radial_var, mask
(4)数据协议代码
脚本位置： pycwr_master/pycwr/io/BaseDataProtocal/AXPT_NRD_Protocol.py
# -*- coding: utf-8 -*-
import numpy as np
SHORT = 'h'
USHORT = 'H'
UCHAR = 'B'
INT = 'i'
UINT = 'I'
FLOAT = 'f'
UCHAR2000 = "B"
class AXPT0164Format(object):
def __init__(self):
self.RadialDataSize = 16104 #径向（Radial）数据大小 == 100 + 8*2000 + 4
self.RadarInfoHeaderBlockPos = 0 #雷达信息（通用+站点）头配置起始
self.TaskConfigurationBlockPos = 32 #任务头配置起始
self.DataBlockPos = 100 #数据起始
self.RadarInfoHeaderBlockSize = 32 #雷达信息头配置大小
self.TaskConfigurationBlockSize = 68 #任务头配置大小
self.flag2Product = {1:'dBT', 2:'dBZ', 3:'V', 4:'W', \
5:'ZDR', 6:'PhiDP', 7:'KDP', 8:'CC'}
self.flag2Product111 = {1:'dBT', 2:'dBZ', 3:'V', 4:'W', 5:'SQI', 6:'CPA', 7:'ZDR', 8:'LDR',\
9:'CC', 10:'PhiDP', 11:'KDP', 12:'CP', 13:'FLAG', 14:'HCL', 15:'CF',\
16:'SNR', 32:'Zc', 33:'Vc', 34:'Wc', 35:'ZDRc'}
self.MomentNum = 2000
self.BaseDataHeader = dict(
RadarInfoHeaderBlock = (
('GenericType', USHORT), # 文件类型 1-基数据文件 2-气象产品文件
('Longitude', SHORT), # 雷达站天线所在位置经度 degree
('Latitude', SHORT), # 雷达站天线所在位置纬度 degree
('BeamWidthVert', UCHAR), # 垂直波束宽度 degree
('BeamWidthHori', UCHAR), # 水平波束宽度 degree
('Height', USHORT), # 天线馈源水平时海拔高度 meter
('SiteName', '22s') # 保留字段(站号)
),
TaskConfigurationBlock = (
('VolumeStartTime', UINT), # 扫描开始秒 seconds 扫描开始时间为UTC标准时间计数,1970年1月1日0时为起始计数基准点
('VolumeStartDate', UINT), # 扫描开始日期 seconds
('NyquistDistance', USHORT), # 理论最大不模糊距离 km
('Azimuth', USHORT), # 方位角 degree
('AzimuthNumber', USHORT), # 方位数
('Elevation', USHORT), # 仰角 degree
('ElevationNumber', USHORT), # 仰角数
('SequenceNumber', USHORT), # 当前方位角内，仰角编号 从0开始
('StartGateRange', USHORT), # 起始库距离 米
('UnitGateLength', USHORT), # 单位库距离 米
('LengthOfGate', USHORT), # 本径向数据块的长度
('Reserved02', '26s'), # 保留字段 （扇区号+系统定正常数+体扫模式）
('RefdBT', USHORT), # 数据回放指针：原始反射率因子
('RefdBZ', USHORT), # 数据回放指针：反射率因子
('RefV', USHORT), # 数据回放指针：径向速度反射率因子
('RefW', USHORT), # 数据回放指针：谱宽
('RefZDR', USHORT), # 数据回放指针：差分反射率因子
('RefPhiDP', USHORT), # 数据回放指针：差分相位
('RefKDP', USHORT), # 数据回放指针：比差分相位
('RefCC', USHORT), # 数据回放指针：相关系数
),
RadialFieldDataBlock = (
('dBT', UCHAR2000), # 数据：原始反射率因子
('dBZ', UCHAR2000), # 数据：反射率因子
('V', UCHAR2000), # 数据：径向速度反射率因子
('W', UCHAR2000), # 数据：谱宽
('ZDR', UCHAR2000), # 数据：差分反射率因子
('PhiDP', UCHAR2000), # 数据：差分相位
('KDP', UCHAR2000), # 数据：比差分相位
('CC', UCHAR2000), # 数据：相关系数
),
)
self.n_variable = len(self.BaseDataHeader['RadialFieldDataBlock'])
dtype_AXPT0164 = AXPT0164Format()
(5)测试主代码和输出
脚本位置：pycwr_master/check.py
pathnow = os.path.split(os.path.realpath(__file__))[0]
import pyximport; pyximport.install()
sys.path.append(f"{pathnow}/")
import pycwr.io.AXPT_NRD as AXPT_NRD
if "__main__" == "__main__":
time0 = time.time()
filename=f"{pathnow}/Z_RADR_I_ZGZ01_20200819000008_O_DOR_DXK_CAR.bin.bz2"
rdata = AXPT_NRD.AXPT01642NRadar( AXPT_NRD.AXPT0164BaseData(filename) )
print(rdata)
for child_var in rdata.fields.keys():
vol_array = rdata.fields[child_var]
print(f' {child_var}: {vol_array.shape}, {np.unique(vol_array)[0:3]} {np.min(vol_array)} ~ {np.max(vol_array)})')
print(f' radar_io: {time.time()-time0} seconds')
读取“Z_RADR_I_ZGZ01_20200819000008_O_DOR_DXK_CAR.bin.bz2”时，
输出：
<pycwr.io.AXPT_NRD.AXPT01642NRadar object at 0x7f6f068c1050>
dBT: (4800, 1383), [-988. -20. -19.58431373] -988.0 ~ 85.58431372549019)
dBZ: (4800, 1383), [-988. -20. -19.58431373] -988.0 ~ 85.58431372549019)
V: (4800, 1383), [-988. -26.23529412 -25.83921569] -988.0 ~ 26.443137254901956)
W: (4800, 1383), [-988. -50. -0.09411765] -988.0 ~ 14.560784313725492)
ZDR: (4800, 1383), [-988. -5. -4.9372549] -988.0 ~ 10.937254901960785)
PhiDP: (4800, 1383), [-988. 0. 0.70980392] -988.0 ~ 179.58039215686276)
KDP: (4800, 1383), [-988. -9.75686275 -9.63529412] -988.0 ~ 20.87843137254902)
CC: (4800, 1383), [-988. 0. 0.02156863] -988.0 ~ 1.0956862745098042)
radar_io: 8.408937215805054 seconds
(6)产品图例
本文参与 腾讯云自媒体分享计划，分享自微信公众号。
原始发表：2021-11-09，如有侵权请联系 cloudcommunity@tencent.com 删除
评论
登录后参与评论
推荐阅读