繼續(xù)雷達(dá)拼圖的工作
本文主要介紹如何將CINRAD雷達(dá)數(shù)據(jù)轉(zhuǎn)換為pyart所需要的格式磕潮。
經(jīng)過查看pyart官方文檔弛矛，可以發(fā)現(xiàn)佣耐，盡管pyart能夠讀取多種雷達(dá)數(shù)據(jù)格式矫夯，但對我國的新一代天氣雷達(dá)基本是無效的，因此需要對pyart進(jìn)行有益的補(bǔ)充，也就是將我們的CINrAD讀取方式加入测蘑。
好了灌危，到此整理一下思路，分以下步驟：
- 查看pyart的格式以及其讀取其他雷達(dá)數(shù)據(jù)的方法
- 查看CINRAD基數(shù)據(jù)的讀取方式
- 綜合前兩步驟碳胳，寫出讀取讀取代碼
查看pyart的格式
隨便打開一個(gè)pyart的例子勇蝙，我們可以發(fā)現(xiàn)pyart將讀取模塊命名為io，我們找一下和雷達(dá)拼圖有關(guān)的代碼
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)
可以發(fā)現(xiàn)在io中有個(gè)read_cfradial函數(shù)挨约，于是順藤摸瓜找到這個(gè)read_cfradial函數(shù)存在于io文件夾下的cfradial.py文件味混。
進(jìn)入下一步，打開文件诫惭，找到read_cfradial函數(shù)翁锡。那就慢慢的看吧，于是又發(fā)現(xiàn)這個(gè)函數(shù)的大部分內(nèi)容都是沒用的夕土，因?yàn)樗拇蟛糠謨?nèi)容都是對雷達(dá)數(shù)據(jù)進(jìn)行整理的馆衔，我們不關(guān)心，因?yàn)镃INRAD的格式和它肯定不同怨绣，因此不需要參考哈踱。
最終看一下這個(gè)函數(shù)返回的結(jié)果吧，見下圖：
我們發(fā)現(xiàn)：
- 返回了一個(gè)
Radar類 - 給類傳遞了很多參數(shù)（暫時(shí)不理會）
到這里就大概有了了解梨熙，這個(gè)讀取的方式也并不是那么的復(fù)雜，不過就是將數(shù)據(jù)讀取進(jìn)來進(jìn)行整理并返回pyart的Radar類刀诬。那如何了解該函數(shù)是如何整理的呢咽扇？
另辟蹊徑的方式就是根據(jù)pyart給出的例子，看看它讀取出來的數(shù)據(jù)具有什么樣的結(jié)構(gòu)不就可以了陕壹？結(jié)論如下：
-
time,range,metadata,scan_type.....,除scan_type是字符串外质欲，其他參數(shù)都是字典 - 每個(gè)參數(shù)的字典結(jié)構(gòu)類似netcdf格式，其中data鍵是數(shù)據(jù)糠馆，其他鍵類似netcdf格式中的attribute嘶伟。
- 數(shù)據(jù)的組成形式，pyart很有意思又碌，它儲存數(shù)據(jù)的方式比較特殊九昧，我們知道雷達(dá)是一條射線一樣發(fā)出輻射的，每次接收回來的數(shù)據(jù)也是一條毕匀。那么它就將這些所有的射線數(shù)據(jù)存儲為二維數(shù)組铸鹰，第一維對應(yīng)于每條射線，第二維對應(yīng)于每條射線上的數(shù)據(jù)皂岔。
當(dāng)然這里說pyart很有意思是在沒有查看CINRAD數(shù)據(jù)格式之前蹋笼。
到這里關(guān)于pyart數(shù)據(jù)結(jié)構(gòu)部分就完事了。下面來看一下CINRAD基數(shù)據(jù)讀取。
CINRAD
寫在最前的參考內(nèi)容
IDL讀取CC雷達(dá)
Python 新一代多普勒天氣雷達(dá)基數(shù)據(jù)可視化
PyCINRAD
CINRAD雷達(dá)產(chǎn)品顯示系統(tǒng)
第一步當(dāng)然是度娘一下CINRAD相關(guān)內(nèi)容剖毯，了解它的數(shù)據(jù)格式等等圾笨，終于讓我發(fā)現(xiàn)了一下可供參考的信息。
終于看到雷達(dá)基數(shù)據(jù)格式了逊谋，我們已經(jīng)知道基數(shù)據(jù)以二進(jìn)制的方式保存擂达，使用python讀取基數(shù)據(jù)是非常方便的，類似于處理Grads數(shù)據(jù)涣狗。直接使用
np.fromfile即可谍婉，當(dāng)然也可以采用如下方式：
with open(file, 'rb') as fopen:
data = np.fromfile(fopen.read(), dtype=***, count=1)
這里的dtype也就是對應(yīng)于基數(shù)據(jù)格式，count則是讀取幾次dtype镀钓，如果為-1則全部讀取穗熬。
對于上圖的基數(shù)據(jù)，我們很容易就可以寫出它的dtype：
header_dtype = np.dtype(
[
('header1', 'u2', 7),
(radar_data', 'u2'),
('header2', 'u2', 6),
...
]
)
當(dāng)然這里也搜索到一些C#的代碼供參考丁溅，可以將其轉(zhuǎn)換為對應(yīng)的python中numpy庫能夠識別的數(shù)據(jù)類型即可唤蔗。
后來進(jìn)一步的搜索發(fā)現(xiàn)已經(jīng)有國內(nèi)的小伙伴完成了部分工作，并且做得非常不錯窟赏，例如PyCINRAD和Python 新一代多普勒天氣雷達(dá)基數(shù)據(jù)可視化妓柜，既然這樣那就可以參考他們的工作來完成相關(guān)的基數(shù)據(jù)讀取部分了。根據(jù)薛志遠(yuǎn)的描述涯穷，雷達(dá)數(shù)據(jù)每個(gè)方位角（射線）的數(shù)據(jù)都包括了文件頭基數(shù)據(jù)部分保留字段棍掐。
讀取的具體代碼就不在這里貼出來了。
最終經(jīng)過整理拷况，將SA/SB型號的雷達(dá)基數(shù)據(jù)轉(zhuǎn)換為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
)
上述代碼沒有經(jīng)過深度整理作煌，比較粗糙，未來經(jīng)過整理也會在Github上開源赚瘦。
下面對結(jié)果進(jìn)行檢驗(yàn)粟誓，通過繪制ppi來檢驗(yàn)。
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')
二者對比相差不大
結(jié)論
完成了SA/SB型號雷達(dá)起意，下面的雷達(dá)型號會比較的快了鹰服。加油！揽咕！
于2019年6月10日23:42:41
同理對于CA/CB雷達(dá)僅僅是雷達(dá)基數(shù)據(jù)格式的不同悲酷，其他處理方式都是相同的。
于2019年6月10日23:59:07
