import numpy as np
from collections import OrderedDict, Counter
from binfile import BinFile
[docs]class NCDCLevel2File(BinFile):
"""
Read a raw Level II file downloaded from NCDC.
Args:
file_name (str): The name of the file to load.
mode (str): Read/write mode of the file. Default value is 'r', the only currently supported option.
"""
_nl2f_header = OrderedDict(
[('data_fmt', '8s'), ('__dummy__', '1s'), ('volume_num', '3s'), ('julian_day', 'i'), ('msecs', 'i'),
('site_id', '4s')]
)
_nl2f_record_meta = OrderedDict(
[('__dummy__', '12s'), ('rec_size', 'h'), ('channel_id', 'b'), ('msg_type', 'b'), ('sequence_id', 'H'),
('msg_julian_date', 'h'), ('msg_msecs', 'i'), ('segment_count', 'h'), ('segment_num', 'h')]
)
_nl2f_record = OrderedDict(
[('id', '4s'), ('data_msecs', 'i'), ('data_julian_day', 'h'), ('radial_num', 'h'), ('azimuth', 'f'),
('compression', 'b'), ('sp', 'b'), ('rlength', 'h'), ('ars', 'b'), ('rs', 'b'), ('elevation_num', 'b'),
('cut', 'b'), ('elevation', 'f'), ('rsbs', 'b'), ('aim', 'b'), ('dcount', 'h'), ('dpb', 'i' * 9)]
)
_nl2f_param = OrderedDict(
[('gate_count', 'h'), ('first_gate', 'h'), ('gate_size', 'h'), ('ref_threshold', 'h'), ('snr_threshold', 'h'),
('__dummy__', 'h'), ('scale', 'f'), ('offset', 'f')]
)
def __init__(self, file_name, mode='r', byteorder='>'):
super(NCDCLevel2File, self).__init__(file_name, mode, byteorder)
self._msg_bytes = 2432
self._cache = {}
if 'r' in mode:
self._read_headers()
self._read_records()
return
def _read_headers(self):
self._header = {}
self._read_block(NCDCLevel2File._nl2f_header, self._header)
return
def _read_records(self):
self._records = []
while not self._ateof():
record = {}
self._read_block(NCDCLevel2File._nl2f_record_meta, record)
if record['msg_type'] == 31:
self._decode_data_message(record)
self._records.append(record)
else:
self._decode_meta_message(record)
self._group_sweeps()
return
def _decode_meta_message(self, record):
# No-op for now
data_start = self._tell()
self._seek(data_start + self._msg_bytes - 28)
return
def _decode_data_message(self, record):
data_start = self._tell()
self._read_block(NCDCLevel2File._nl2f_record, record)
self._seek(data_start + record['dpb'][0] + 40)
record['vcp'] = self._read('h')
self._read_data(data_start, record['dpb'], record)
self._seek(data_start + (record['rec_size'] - 8) * 2)
return
def _read_data(self, data_start, offsets, record, debug=False):
for dp in range(3, len(offsets)):
if offsets[dp] <= 0:
continue
param = {}
self._seek(data_start + offsets[dp] + 1)
param_name = self._read('7s')
self._read_block(NCDCLevel2File._nl2f_param, param)
dtype = 'H' if param_name in ['PHI'] else 'B'
param['data'] = np.array(self._read("%d%s" % (param['gate_count'], dtype)))
param['data'] = (param['data'] - param['offset']) / param['scale']
record[param_name] = param
return
def _group_sweeps(self):
def _circle_diff(az1, az2):
if az1 > 270 and az2 < 90:
az1 -= 360
elif az1 < 90 and az2 > 270:
az2 -= 360
return az1 - az2
sweep_breaks = [0]
sweep_breaks.extend(idx + 1 for idx in range(len(self._records) - 1)
if _circle_diff(self._records[idx + 1]['azimuth'], self._records[idx]['azimuth']) > 2)
sweep_breaks.append(None)
self._sweeps = {}
self._ranges = {}
self._azimuths = {}
self._elevations = {}
for idx in range(len(sweep_breaks) - 1):
sw_start = sweep_breaks[idx]
sw_end = sweep_breaks[idx + 1]
azimuths = np.array([r['azimuth'] for r in self._records[sw_start:sw_end]])
n_azim = azimuths.shape[0]
if n_azim > 540:
azimuths = np.floor(azimuths * 2) / 2
else:
azimuths = np.floor(azimuths)
elevation = Counter(r['elevation'] for r in self._records[sw_start:sw_end]).most_common(1)[0][0]
# self._elevations.append(elevation)
params = [k for k, v in self._records[sw_start].iteritems() if type(v) in [dict]]
for param in params:
n_range = self._records[sw_start][param]['gate_count']
r_first = self._records[sw_start][param]['first_gate']
dr_gate = self._records[sw_start][param]['gate_size']
ranges = np.arange(n_range) * dr_gate + r_first
swp_array = np.empty((n_azim, n_range))
try:
self._sweeps[param].append(swp_array)
self._azimuths[param].append(azimuths)
self._ranges[param].append(ranges)
self._elevations[param].append(elevation)
except KeyError:
self._sweeps[param] = [swp_array]
self._azimuths[param] = [azimuths]
self._ranges[param] = [ranges]
self._elevations[param] = [elevation]
for jdy, rad in enumerate(self._records[sw_start:sw_end]):
swp_array[jdy] = rad[param]['data']
for param, elevs in self._elevations.iteritems():
self._elevations[param] = np.round(elevs, 1)
return
[docs] def __getitem__(self, var_name):
"""
Retrieve data from the file.
Args:
var_name (str): Name of the variable to retrieve from the file. Acceptable values are 'REF' for
reflectivity and 'VEL' for radial velocity.
Returns:
A three-dimensional numpy array (NTILT :math:`\\times` NAZIM :math:`\\times` NRANGE)
"""
try:
return self._sweeps[var_name]
except KeyError:
raise ValueError("Variable '%s' not found in file." % var_name)
[docs] def get_elevations(self, var_name):
"""
Get the elevation angles for a particular variable.
Args:
var_name (str): The variable name for which to retrieve the elevation angles.
Returns:
A 1-dimensional numpy array containing the elevation angles in degrees.
"""
try:
return self._elevations[var_name]
except KeyError:
raise ValueError("Variable '%s' not found in file." % var_name)
[docs] def get_coords(self, var_name, sweep_no):
"""
Get the azimuth and range for a sweep in the data file.
Args:
var_name (str): The variable name for which to retrieve the azimuth and range.
sweep_no (int): The tilt for which to retrieve azimuth and range.
Returns:
A tuple of 1-dimensional numpy arrays of azimuth (in degrees) and range (in meters).
"""
return self._azimuths[var_name][sweep_no], self._ranges[var_name][sweep_no]
if __name__ == "__main__":
name = "/data6/tsupinie/20110524/wsr-88d/raw/KGLD20110524_212909_V03"
nl2f = NCDCLevel2File(name)
ref = nl2f['REF'][1]
azi_ref, rng_ref = nl2f.get_coords('REF', 1)
azi_ref = np.append(azi_ref, azi_ref[0])
vel = nl2f['VEL'][0]
azi_vel, rng_vel = nl2f.get_coords('VEL', 0)
azi_vel = np.append(azi_vel, azi_vel[0])