forked from YvZheng/pycwr
-
Notifications
You must be signed in to change notification settings - Fork 0
/
SCFile.py
414 lines (366 loc) · 17 KB
/
SCFile.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
# -*- coding: utf-8 -*-
import numpy as np
from .BaseDataProtocol.SCProtocol import dtype_sc
from .util import _prepare_for_read, _unpack_from_buf, make_time_unit_str, get_radar_sitename
import pandas as pd
import datetime
from ..core.NRadar import NuistRadar
from ..configure.pyart_config import get_metadata, get_fillvalue
from ..configure.default_config import CINRAD_field_mapping, _LIGHT_SPEED
from ..core.PyartRadar import Radar
from netCDF4 import date2num
class SCBaseData(object):
"""
解码SC/CD 1.0的数据格式
"""
def __init__(self, filename):
super(SCBaseData, self).__init__()
self.filename = filename
self.fid = _prepare_for_read(self.filename) ##判断是否是压缩文件
buf_header = self.fid.read(dtype_sc.BaseDataHeaderSize) ##取出header的buf
self.header = self._parse_BaseDataHeader(buf_header)
self.MaxV = self.header['LayerParam']['MaxV'][0]/100. ##??可能会存在问题,如果不同仰角采用不用的PRF
self._check_sc_basedata()
self.fid.seek(dtype_sc.BaseDataHeaderSize, 0) ##移动到径向数据的位置
self.radial = self._parse_radial()
self.fid.close()
def _check_sc_basedata(self):
"""检查雷达数据是否完整"""
buf_radial_data = self.fid.read()
assert len(buf_radial_data) == self.nrays * dtype_sc.PerRadialSize, "SC basedata size has problems!"
return
def _parse_BaseDataHeader(self, buf_header):
"""
:param buf_header: 只包含头文件的buf
:return:
"""
BaseDataHeader_dict = {}
##解码雷达站点信息
BaseDataHeader_dict['RadarSite'], _ = _unpack_from_buf(buf_header,\
dtype_sc.RadarSitePos,dtype_sc.BaseDataHeader['RadarSite'])
##解码雷达性能参数
BaseDataHeader_dict['RadarPerformanceParam'], _ = _unpack_from_buf(buf_header,\
dtype_sc.RadarPerformanceParamPos,dtype_sc.BaseDataHeader['RadarPerformanceParam'])
##解码观测参数
BaseDataHeader_dict['RadarObserationParam_1'], _ = _unpack_from_buf(buf_header, \
dtype_sc.RadarObserationParamPos_1, dtype_sc.BaseDataHeader['RadarObserationParam_1'])
##目前仅支持VOL扫描格式
assert BaseDataHeader_dict['RadarObserationParam_1']['stype'] > 100, "only vol support!"
self.nsweeps = BaseDataHeader_dict['RadarObserationParam_1']['stype'] - 100
##解码不同仰角的观测参数
BaseDataHeader_dict['LayerParam'] = np.frombuffer(buf_header, \
dtype_sc.BaseDataHeader['LayerParamX30'],count=self.nsweeps, offset=dtype_sc.LayerParamPos)
##解码其余一些观测参数
BaseDataHeader_dict['RadarObserationParam_2'], _ = _unpack_from_buf(buf_header,\
dtype_sc.RadarObserationParamPos_2, dtype_sc.BaseDataHeader['RadarObserationParam_2'])
self.nrays = np.sum(BaseDataHeader_dict['LayerParam']['recordnumber'])
self.sweep_end_ray_index_add1 = np.cumsum(BaseDataHeader_dict['LayerParam']['recordnumber']) ##python格式的结束
self.sweep_start_ray_index = self.sweep_end_ray_index_add1 - BaseDataHeader_dict['LayerParam']['recordnumber']
return BaseDataHeader_dict
def _parse_radial(self):
radial = []
for isweep in range(self.nsweeps):
MaxV = self.header['LayerParam']['MaxV'][isweep] / 100.
for _ in range(self.header['LayerParam']['recordnumber'][isweep]):
buf_radial = self.fid.read(dtype_sc.PerRadialSize)
radial.append(self._parse_radial_single(buf_radial, 0, MaxV, -1))
return radial
def _parse_radial_single(self, buf_radial, start_pos, MaxV, num_bins=-1):
"""
:param buf_radial:
:param start_pos: 开始的pos
:param num_bins: 库数
:return:
"""
radial_dict = {}
radial_dict_tmp, size_tmp = _unpack_from_buf(buf_radial, start_pos, dtype_sc.RadialHeader())
radial_dict.update(radial_dict_tmp)
radial_dict['fields'] = {}
RadialData = np.frombuffer(buf_radial, dtype_sc.RadialData(),\
count=num_bins, offset=start_pos+size_tmp)
radial_dict['fields']['dBZ'] = np.where(RadialData['dBZ'] != 0,\
(RadialData['dBZ'].astype(int) - 64)/2., np.nan).astype(np.float32)
radial_dict['fields']['dBT'] = np.where(RadialData['dBT'] != 0, \
(RadialData['dBT'].astype(int) - 64) / 2., np.nan).astype(np.float32)
radial_dict['fields']['V'] = np.where(RadialData['V'] != 0, \
MaxV * (RadialData['V'].astype(int) - 128) / 128., np.nan).astype(np.float32)
radial_dict['fields']['W'] = np.where(RadialData['W'] != 0, \
MaxV * RadialData['W'].astype(int) /256., np.nan).astype(np.float32)
return radial_dict
def get_nyquist_velocity(self):
"""get nyquist vel per ray
获取每根径向的不模糊速度
:return:(nRays)
"""
nyquist_velocity = np.concatenate([np.array([self.header['LayerParam']['MaxV'][isweep] / \
100.] * self.header['LayerParam']['recordnumber'][isweep]) for \
isweep in range(self.nsweeps)])
return nyquist_velocity.astype(np.float32)
def get_unambiguous_range(self):
"""
获取每根径向的不模糊距离
:return:(nRays)
"""
return np.concatenate([np.array([self.header['LayerParam']['MaxL'][isweep] *10 \
] * self.header['LayerParam']['recordnumber'][isweep]) for \
isweep in range(self.nsweeps)])
def get_scan_time(self):
"""
获取每根径向的扫描时间
:return:(nRays)
"""
Start_params = self.header['RadarObserationParam_1']
End_params = self.header['RadarObserationParam_2']
start_time = datetime.datetime(year=Start_params['syear'], month=Start_params['smonth'],
day=Start_params['sday'], hour=Start_params['shour'],
minute=Start_params['sminute'], second=Start_params['ssecond'])
end_time = datetime.datetime(year=End_params['Eyear'], month=End_params['Emonth'],
day=End_params['Eday'], hour=End_params['Ehour'],
minute=End_params['Eminute'], second=End_params['Esecond'])
return pd.date_range(start_time, end_time, periods=self.nrays).to_pydatetime() - datetime.timedelta(hours=8)
def get_sweep_end_ray_index(self):
"""
获取每个sweep的结束的index,包含在内
:return:(nsweep)
"""
return self.sweep_end_ray_index_add1 - 1
def get_sweep_start_ray_index(self):
"""
获取每个sweep的开始的index
:return:(nsweep)
"""
return self.sweep_start_ray_index
def get_rays_per_sweep(self):
"""
获取每个sweep的径向数
:return:(nsweep)
"""
return self.sweep_end_ray_index_add1 - self.sweep_start_ray_index
def get_azimuth(self):
"""
获取每根径向的方位角
:return:(nRays)
"""
return np.array([(iray['sStrAz'] + iray['sEndAz'])*180./65536 for iray in self.radial])
def get_elevation(self):
"""
获取每根径向的仰角
:return: (nRays)
"""
return np.array([(iray['sStrEl'] + iray['sEndEl'])*180./65536 for iray in self.radial])
def get_latitude_longitude_altitude_frequency(self):
"""
获取经纬度高度,雷达频率
:return:lat, lon, alt, frequency
"""
return self.header['RadarSite']['latitudevalue']/100., self.header['RadarSite']['longitudevalue']/100., \
self.header['RadarSite']['height'] / 1000., 2.765
def get_scan_type(self):
"""
获取扫描的类型
:return:
"""
if self.header['RadarObserationParam_1']['stype'] == 1:
return "rhi"
else:
return "ppi"
def get_sitename(self):
return get_radar_sitename(self.filename)
class SC2NRadar(object):
"""到NusitRadar object 的桥梁"""
def __init__(self, SC):
self.SC = SC
self.radial = self.SC.radial
self.azimuth = self.get_azimuth()
self.elevation = self.get_elevation()
self.sweep_start_ray_index = self.get_sweep_start_ray_index()
self.sweep_end_ray_index = self.get_sweep_end_ray_index()
self.nrays = self.SC.nrays
self.nsweeps = self.SC.nsweeps
self.scan_type = self.SC.get_scan_type()
self.latitude, self.longitude, self.altitude, self.frequency = \
self.SC.get_latitude_longitude_altitude_frequency()
self.bins_per_sweep = self.get_nbins_per_sweep()
self.max_bins = self.bins_per_sweep.max()
self.range = self.get_range_per_radial(self.max_bins)
self.fields = self._get_fields()
self.sitename = self.SC.get_sitename()
def get_azimuth(self):
"""
获取每根径向的方位角
:return:(nRays)
"""
return self.SC.get_azimuth()
def get_elevation(self):
"""
获取每根径向的仰角
:return: (nRays)
"""
return self.SC.get_elevation()
def get_rays_per_sweep(self):
"""
获取每个sweep的径向数
:return:(nsweep)
"""
return (self.SC.header['LayerParam']['recordnumber']).astype(int)
def get_scan_time(self):
"""
获取每根径向的扫描时间
:return:(nRays)
"""
return self.SC.get_scan_time()
def get_nyquist_velocity(self):
"""get nyquist vel per ray
获取每根径向的不模糊速度
:return:(nRays)
"""
return self.SC.get_nyquist_velocity()
def get_unambiguous_range(self):
"""
获取每根径向的不模糊距离
:return:(nRays)
"""
return self.SC.get_unambiguous_range()
def get_sweep_end_ray_index(self):
"""
获取每个sweep的结束的index,包含在内
:return:(nsweep)
"""
return self.SC.sweep_end_ray_index_add1 - 1
def get_sweep_start_ray_index(self):
"""
获取每个sweep的开始的index
:return:(nsweep)
"""
return self.SC.sweep_start_ray_index
def get_nbins_per_sweep(self):
"""
确定每个sweep V探测的库数
:return:
"""
return np.array([self.radial[idx]['fields']['dBZ'].size for idx in self.sweep_start_ray_index])
def get_range_per_radial(self, length):
"""
确定径向每个库的距离
:param length:
:return:
"""
Resolution = self.SC.header['LayerParam']['binWidth'][0]/10.
return np.linspace(Resolution, Resolution * length, length)
def _get_fields(self):
"""将所有的field的数据提取出来"""
fields = {}
field_keys = self.radial[0]['fields'].keys()
for ikey in field_keys:
fields[ikey] = np.array([(iray['fields'][ikey]).ravel() for iray in self.radial])
return fields
def get_NRadar_nyquist_speed(self):
"""array shape (nsweeps)"""
return self.SC.header['LayerParam']['MaxV'] / 100.
def get_NRadar_unambiguous_range(self):
"""array shape (nsweeps)"""
return self.SC.header['LayerParam']['MaxL'] * 10.
def get_fixed_angle(self):
return self.SC.header['LayerParam']['Swangles'] / 100.
def ToNuistRadar(self):
"""将WSR98D数据转为Nuist Radar的数据格式"""
return NuistRadar(fields=self.fields, scan_type=self.scan_type, time=self.get_scan_time(), \
range=self.range, azimuth=self.azimuth, elevation=self.elevation, latitude=self.latitude, \
longitude=self.longitude, altitude=self.altitude,
sweep_start_ray_index=self.sweep_start_ray_index, \
sweep_end_ray_index=self.sweep_end_ray_index, fixed_angle=self.get_fixed_angle(), \
bins_per_sweep=self.bins_per_sweep, nyquist_velocity=self.get_NRadar_nyquist_speed(), \
frequency=self.frequency, unambiguous_range=self.get_NRadar_unambiguous_range(), \
nrays=self.nrays, nsweeps=self.nsweeps, sitename = self.sitename)
def ToPyartRadar(self):
dts = self.get_scan_time()
units = make_time_unit_str(min(dts))
time = get_metadata('time')
time['units'] = units
time['data'] = date2num(dts, units).astype('float32')
# range
_range = get_metadata('range')
# assume that the number of gates and spacing from the first ray is
# representative of the entire volume
_range['data'] = self.range
_range['meters_to_center_of_first_gate'] = self.range[0]
_range['meters_between_gates'] = self.range[0]
latitude = get_metadata('latitude')
longitude = get_metadata('longitude')
altitude = get_metadata('altitude')
latitude['data'] = np.array([self.latitude], dtype='float64')
longitude['data'] = np.array([self.longitude], dtype='float64')
altitude['data'] = np.array([self.altitude], dtype='float64')
metadata = get_metadata('metadata')
metadata['original_container'] = 'CINRAD/SC'
metadata['site_name'] = self.sitename
metadata['radar_name'] = "CINRAD/SC"
sweep_start_ray_index = get_metadata('sweep_start_ray_index')
sweep_end_ray_index = get_metadata('sweep_end_ray_index')
sweep_start_ray_index['data'] = self.sweep_start_ray_index
sweep_end_ray_index['data'] = self.sweep_end_ray_index
sweep_number = get_metadata('sweep_number')
sweep_number['data'] = np.arange(self.nsweeps, dtype='int32')
scan_type = self.scan_type
sweep_mode = get_metadata('sweep_mode')
if self.scan_type == "ppi":
sweep_mode['data'] = np.array(self.nsweeps * ['azimuth_surveillance'], dtype='S')
elif self.scan_type == "rhi":
sweep_mode['data'] = np.array(self.nsweeps * ['rhi'], dtype='S')
else:
sweep_mode['data'] = np.array(self.nsweeps * ['sector'], dtype='S')
# elevation
elevation = get_metadata('elevation')
elevation['data'] = self.elevation
# azimuth
azimuth = get_metadata('azimuth')
azimuth['data'] = self.azimuth
# fixed_angle
fixed_angle = get_metadata('fixed_angle')
fixed_angle['data'] = self.get_fixed_angle()
# instrument_parameters
instrument_parameters = self._get_instrument_parameters()
# fields
fields = {}
for field_name_abbr in self.fields.keys():
field_name = CINRAD_field_mapping[field_name_abbr]
if field_name is None:
continue
field_dic = get_metadata(field_name)
field_dic['data'] = np.ma.masked_array(self.fields[field_name_abbr],\
mask=np.isnan(self.fields[field_name_abbr]), fill_value=get_fillvalue())
field_dic['_FillValue'] = get_fillvalue()
fields[field_name] = field_dic
return Radar(time, _range, fields, metadata, scan_type,
latitude, longitude, altitude,
sweep_number, sweep_mode, fixed_angle, sweep_start_ray_index,
sweep_end_ray_index,
azimuth, elevation,
instrument_parameters=instrument_parameters)
def _get_instrument_parameters(self):
""" Return a dictionary containing instrument parameters. """
# pulse width
pulse_width = get_metadata('pulse_width')
pulse_width['data'] = self.range[0] / _LIGHT_SPEED # nanosec->sec
# assume that the parameters in the first ray represent the beam widths,
# bandwidth and frequency in the entire volume
wavelength_hz = self.frequency * 10 ** 9
# radar_beam_width_h
radar_beam_width_h = get_metadata('radar_beam_width_h')
radar_beam_width_h['data'] = np.array([1., ], dtype='float32')
# radar_beam_width_v
radar_beam_width_v = get_metadata('radar_beam_width_v')
radar_beam_width_v['data'] = np.array([1., ], dtype='float32')
# frequency
frequency = get_metadata('frequency')
frequency['data'] = np.array([wavelength_hz, ], dtype='float32')
instrument_parameters = {
'pulse_width': pulse_width,
'radar_beam_width_h': radar_beam_width_h,
'radar_beam_width_v': radar_beam_width_v,
'frequency': frequency, }
# nyquist velocity if defined
nyquist_velocity = get_metadata('nyquist_velocity')
nyquist_velocity['data'] = self.get_nyquist_velocity()
instrument_parameters['nyquist_velocity'] = nyquist_velocity
return instrument_parameters