CINRAD数据处理

继续雷达拼图的工作
本文主要介绍如何将CINRAD雷达数据转换为pyart所需要的格式。

经过查看pyart官方文档，可以发现，尽管pyart能够读取多种雷达数据格式，但对我国的新一代天气雷达基本是无效的，因此需要对pyart进行有益的补充，也就是将我们的CINrAD读取方式加入。
好了，到此整理一下思路，分以下步骤：

查看pyart的格式以及其读取其他雷达数据的方法
查看CINRAD基数据的读取方式
综合前两步骤，写出读取读取代码

查看pyart的格式

随便打开一个pyart的例子，我们可以发现pyart将读取模块命名为io，我们找一下和雷达拼图有关的代码

XSAPR_SW_FILE = 'swx_20120520_0641.nc'
XSAPR_SE_FILE = 'sex_20120520_0641.nc'
radar_sw = pyart.io.read_cfradial(XSAPR_SW_FILE)
radar_se = pyart.io.read_cfradial(XSAPR_SE_FILE)

可以发现在io中有个read_cfradial函数，于是顺藤摸瓜找到这个read_cfradial函数存在于io文件夹下的cfradial.py文件。

进入下一步，打开文件，找到read_cfradial函数。那就慢慢的看吧，于是又发现这个函数的大部分内容都是没用的，因为它的大部分内容都是对雷达数据进行整理的，我们不关心，因为CINRAD的格式和它肯定不同，因此不需要参考。

最终看一下这个函数返回的结果吧，见下图：

image.png

我们发现：

返回了一个Radar类
给类传递了很多参数（暂时不理会）

到这里就大概有了了解，这个读取的方式也并不是那么的复杂，不过就是将数据读取进来进行整理并返回pyart的Radar类。那如何了解该函数是如何整理的呢？

另辟蹊径的方式就是根据pyart给出的例子，看看它读取出来的数据具有什么样的结构不就可以了？结论如下：

time,range,metadata,scan_type.....,除scan_type是字符串外，其他参数都是字典
每个参数的字典结构类似netcdf格式，其中data键是数据，其他键类似netcdf格式中的attribute。
数据的组成形式，pyart很有意思，它储存数据的方式比较特殊，我们知道雷达是一条射线一样发出辐射的，每次接收回来的数据也是一条。那么它就将这些所有的射线数据存储为二维数组，第一维对应于每条射线，第二维对应于每条射线上的数据。

当然这里说pyart很有意思是在没有查看CINRAD数据格式之前。

到这里关于pyart数据结构部分就完事了。下面来看一下CINRAD基数据读取。

CINRAD

写在最前的参考内容
IDL读取CC雷达
 Python 新一代多普勒天气雷达基数据可视化
 PyCINRAD
CINRAD雷达产品显示系统

第一步当然是度娘一下CINRAD相关内容，了解它的数据格式等等，终于让我发现了一下可供参考的信息。

image.png

终于看到雷达基数据格式了，我们已经知道基数据以二进制的方式保存，使用python读取基数据是非常方便的，类似于处理Grads数据。直接使用np.fromfile即可，当然也可以采用如下方式：

with open(file, 'rb') as fopen:
    data = np.fromfile(fopen.read(), dtype=***, count=1)

这里的dtype也就是对应于基数据格式，count则是读取几次dtype，如果为-1则全部读取。
对于上图的基数据，我们很容易就可以写出它的dtype：

header_dtype = np.dtype(
    [
    ('header1', 'u2', 7),
    (radar_data', 'u2'),
    ('header2', 'u2', 6),
    ...
    ]
    )

当然这里也搜索到一些C#的代码供参考，可以将其转换为对应的python中numpy库能够识别的数据类型即可。

后来进一步的搜索发现已经有国内的小伙伴完成了部分工作，并且做得非常不错，例如PyCINRAD和Python 新一代多普勒天气雷达基数据可视化，既然这样那就可以参考他们的工作来完成相关的基数据读取部分了。根据薛志远的描述，雷达数据每个方位角（射线）的数据都包括了文件头基数据部分保留字段。

image.png

读取的具体代码就不在这里贴出来了。
最终经过整理，将SA/SB型号的雷达基数据转换为pyart的Radar类的代码如下：

import os
import re
import datetime
import pyart
from pyart.core import Radar


def __check_base_type(self, filein):
        # Z9517_BASE_SA_20190531_082100.bin
        # D:\Download\radar\Z_RADR_I_Z9010_20190521010000_O_DOR_SA_CAP.bin.bz2
        type_re_math = '_SA_|_SB_|_SC_|_CA_|_CB_|_CC_|_CD_'
        station_re_math = '_Z\d{4}_|Z\d{4}_|_Z\d{4}'
        time_re_math = '_\d{14}_|_\d{8}_|_\d{6}'
        filename = os.path.basename(filein)
        tmp_search = re.findall(type_re_math, filename)
        if tmp_search and len(tmp_search) == 1:
            radar_base_data_type = tmp_search[0].replace('_', '')
        tmp_search = re.findall(station_re_math, filename)
        if tmp_search and len(tmp_search) == 1:
            radar_station_id = tmp_search[0].replace('_', '')
        tmp_search = re.findall(time_re_math, filename)
        if tmp_search:
            if len(tmp_search) == 1:
                radar_scan_time_str = tmp_search[0].replace('_', '')
                radar_scan_time = datetime.datetime.strptime(
                    radar_scan_time_str, '%Y%m%d%H%M%S'
                    ) + datetime.timedelta(hours=8)
            elif len(tmp_search) == 2:
                radar_scan_time_str = ''.join([i.replace('_', '') for i in tmp_search])
                radar_scan_time = datetime.datetime.strptime(
                    radar_scan_time_str, '%Y%m%d%H%M%S'
                    ) 
        return radar_station_id, radar_base_data_type, radar_scan_time

def __prepare_work(self, filein):
        radar_station_id, radar_base_data_type, radar_scan_time = self.__check_base_type(filein)
        radar_station_info = self.stainfo.get(radar_station_id)
        if radar_station_info is not None:
            radar_station_name = radar_station_info[0]
            radar_longitude = radar_station_info[1]
            radar_latitude = radar_station_info[2]
            radar_altitude = radar_station_info[-1]
        else:
            radar_station_name = radar_station_info[0]
            radar_longitude = radar_station_info[1]
            radar_latitude = radar_station_info[2]
            radar_altitude = radar_station_info[-1]

        data_prepare = {
            'time': dict(
                comment='time ozone is BEIJING!',
                calendar = 'standard',
                units = 'seconds since 1970-01-01 00:00',
                standard_name='time',
                long_name='time in seconds since volume start',
                data=[date2num(radar_scan_time,'seconds since 1970-01-01 00:00','standard')]),
            'radar_lon' : dict(
                units = 'degrees_east',
                standard_name = 'Longitude',
                data=np.ma.array([radar_longitude], mask=0)
                ),
            'radar_lat':dict(
                units = 'degrees_north',
                standard_name = 'Latitude',
                data = np.ma.array([radar_latitude], mask=0)
                ),
            'radar_alt':dict(
                units='meters',
                standard_name = 'Altitude',
                data=np.ma.array([radar_altitude], mask=0)
                ),
            'metadata':dict(
                radar_station_name = radar_station_name,
                radar_station_id = radar_station_id,
                radar_base_data_type = radar_base_data_type,
                radar_scan_time = radar_scan_time,
                radar_longitude = radar_longitude,
                radar_latitude = radar_latitude,
                radar_altitude = radar_altitude,
                )
              }
        return data_prepare

def __radar_SA_SB(self, filein):
        data_prepare = self.__prepare_work(filein)
        with self.__return_file_object(filein) as fopen:
            data = np.frombuffer(fopen.read(), dtype=self.SAB_dtype)
        scan_type = 'ppi'
        data_prepare['time']['data'] = data_prepare['time']['data'] * data['azimuth'].size
        azimuth = dict(
            data = (data['azimuth'] / 8.) * (180/4096),
            units = 'degrees',
            standard_name='beam_azimuth_angle',
            long_name='azimuth_angle_from_true_north'
            )
        elevation = dict(
            units='degrees',
            comment='Elevation of antenna relative to the horizontal plane',
            standard_name='beam_elevation_angle',
            long_name='elevation_angle_from_horizontal_plane',
            data=(data['elevation'] / 8.) * (180/4096)
            )

        metadata_new = dict(
            nyquist_vel = data['nyquist_vel'][0]/100,
            instrument_name='China CINRAD'
            )
        data_prepare['metadata'].update(metadata_new)

        _range=dict(
            units='meters',
            meters_between_gates=50,
            long_name='range_to_measurement_volume',
            standard_name='projection_range_coordinate',
            meters_to_center_of_first_gate= 0,
            data = np.ma.array(np.arange(0,460,1), mask=False)
            )

        sweep_start_ray_index = dict(
            data=np.ma.array(np.where(np.logical_or(data['radial_state']==3, data['radial_state']==0))[0], mask=False),
            units='count',
            long_name='index of first ray in sweep, 0-based',
            )
        sweep_end_ray_index = dict(
            data=np.ma.array(np.where(np.logical_or(data['radial_state']==4, data['radial_state']==2))[0], mask=False),
            units='count',
            long_name='index of last ray in sweep, 0-based',
            )

        fix_angle = dict(
            units='degrees',
            long_name='target_angle_for_sweep',
            standard_name = 'target_fixed_angle',
            data=np.ma.array(np.unique(elevation['data']), mask=False)
            )

        sweep_number = dict(
            units='count', 
            long_name='sweep_number',
            data=np.ma.array(np.arange(len(sweep_start_ray_index['data'])), mask=False)
            )

        sweep_mode = dict(
            units='uniteless',
            long_name='sweep_mode',
            comment= 'Options are:"sector","coplane",rhi","vertical_pointing","idle","azimuth_surveillance","elevation_surveillance","sunscan","pointing","manual_ppi","manual_rhi"',
            data = np.ma.array(['manual_ppi'] * len(sweep_start_ray_index['data']), mask=False)
            )

        fields = dict(
            reflectivity_horizontal=dict(
                _FillValue=0,
                units='dBZ',
                long_name='equivalent_reflectivity_factor',
                standard_name='equivalent_reflectivity_factor',
                valid_max=90,
                valid_min=0,
                data=(np.ma.masked_where(data['r']==0, data['r']) -2 ) /2 -32
                )
            )


        return Radar(
            data_prepare['time'], _range,
            fields, data_prepare['metadata'], scan_type,
            data_prepare['radar_lat'],  data_prepare['radar_lon'], data_prepare['radar_alt'],
            sweep_number, sweep_mode, fix_angle, sweep_start_ray_index, sweep_end_ray_index,
            azimuth, elevation
            )

上述代码没有经过深度整理，比较粗糙，未来经过整理也会在Github上开源。
下面对结果进行检验，通过绘制ppi来检验。

fig = plt.figure(figsize=[8, 8])
display = pyart.graph.RadarMapDisplay(radar)
display.plot_ppi_map('reflectivity_horizontal', sweep=0, resolution='50m',
                         vmin=-8, vmax=64, 
                         projection=ccrs.PlateCarree(), cmap='pyart_LangRainbow12')
plt.savefig('../test.png')

test.png

image.png

二者对比相差不大

结论

完成了SA/SB型号雷达，下面的雷达型号会比较的快了。加油！！
于2019年6月10日23:42:41

同理对于CA/CB雷达仅仅是雷达基数据格式的不同，其他处理方式都是相同的。
于2019年6月10日23:59:07