"""
Utilities for reading gamic hdf5 files.
"""
# TODO to move out of aux_io namespace:
# * unit tests - need small sample files.
# * auto-detect file type with pyart.io.read function
# * move to pyart.io namespace
import datetime
import warnings
import numpy as np
from ..config import FileMetadata, get_fillvalue
from ..io.common import make_time_unit_str, _test_arguments
from ..core.radar import Radar
try:
from .gamicfile import GAMICFile
_H5PY_AVAILABLE = True
except ImportError:
_H5PY_AVAILABLE = False
from ..exceptions import MissingOptionalDependency
LIGHT_SPEED = 2.99792458e8 # speed of light in meters per second
[docs]def read_gamic(filename, field_names=None, additional_metadata=None,
file_field_names=False, exclude_fields=None,
include_fields=None, valid_range_from_file=True,
units_from_file=True, pulse_width=None, **kwargs):
"""
Read a GAMIC hdf5 file.
Parameters
----------
filename : str
Name of GAMIC HDF5 file to read data from.
field_names : dict, optional
Dictionary mapping field names in the file names to radar field names.
Unlike other read functions, fields not in this dictionary or having a
value of None are still included in the radar.fields dictionary, to
exclude them use the `exclude_fields` parameter. Fields which are
mapped by this dictionary will be renamed from key to value.
additional_metadata : dict of dicts, optional
This parameter is not used, it is included for uniformity.
file_field_names : bool, optional
True to force the use of the field names from the file in which
case the `field_names` parameter is ignored. False will use to
`field_names` parameter to rename fields.
exclude_fields : list or None, optional
List of fields to exclude from the radar object. This is applied
after the `file_field_names` and `field_names` parameters. Set
to None to include all fields specified by include_fields.
include_fields : list or None, optional
List of fields to include from the radar object. This is applied
after the `file_field_names` and `field_names` parameters. Set
to None to include all fields not specified by exclude_fields.
valid_range_from_file : bool, optional
True to extract valid range (valid_min and valid_max) for all
field from the file when they are present. False will not extract
these parameters.
units_from_file : bool, optional
True to extract the units for all fields from the file when available.
False will not extract units using the default units for the fields.
pulse_width : list or None,
Mandatory for gamic radar processors which have pulsewidth enums.
pulse_width should contain the pulsewidth' in us.
Returns
-------
radar : Radar
Radar object.
"""
# check that h5py is available
if not _H5PY_AVAILABLE:
raise MissingOptionalDependency(
"h5py is required to use read_gamic but is not installed")
# test for non empty kwargs
_test_arguments(kwargs)
# create metadata retrieval object
filemetadata = FileMetadata('gamic', field_names, additional_metadata,
file_field_names, exclude_fields,
include_fields)
# Open HDF5 file and get handle
gfile = GAMICFile(filename)
# verify that all scans are present in file
assert gfile.is_file_complete()
# latitude, longitude and altitude
latitude = filemetadata('latitude')
longitude = filemetadata('longitude')
altitude = filemetadata('altitude')
latitude['data'] = gfile.where_attr('lat', 'float64')
longitude['data'] = gfile.where_attr('lon', 'float64')
altitude['data'] = gfile.where_attr('height', 'float64')
# metadata
metadata = filemetadata('metadata')
metadata['original_container'] = 'GAMIC-HDF5'
what_mapping = {
'version': 'gamic_version', 'sets_scheduled': 'sets_scheduled',
'object': 'gamic_object', 'date': 'gamic_date'}
for gamic_key, metadata_key in what_mapping.items():
if gfile.is_attr_in_group('what', gamic_key):
metadata[metadata_key] = gfile.raw_group_attr('what', gamic_key)
how_keys = ['software', 'template_name', 'site_name', 'host_name',
'azimuth_beam', 'elevation_beam', 'sdp_name', 'sw_version',
'sdp_version', 'simulated']
for key in how_keys:
if gfile.is_attr_in_group('how', key):
metadata[key] = gfile.raw_group_attr('how', key)
# sweep_start_ray_index, sweep_end_ray_index
sweep_start_ray_index = filemetadata('sweep_start_ray_index')
sweep_end_ray_index = filemetadata('sweep_end_ray_index')
sweep_start_ray_index['data'] = gfile.start_ray.astype('int32')
sweep_end_ray_index['data'] = gfile.end_ray.astype('int32')
# sweep number
sweep_number = filemetadata('sweep_number')
try:
sweep_number['data'] = gfile.what_attrs('set_idx', 'int32')
except KeyError:
sweep_number['data'] = np.arange(gfile.nsweeps, dtype='int32')
# sweep_type
scan_type = gfile.raw_scan0_group_attr('what', 'scan_type').lower()
if hasattr(scan_type, 'decode'):
scan_type = scan_type.decode('utf-8')
# check that all scans in the volume are the same type
if not gfile.is_file_single_scan_type():
raise NotImplementedError('Mixed scan_type volume.')
if scan_type not in ['ppi', 'rhi']:
message = "Unknown scan type: %s, reading as RHI scans." % (scan_type)
warnings.warn(message)
scan_type = 'rhi'
# sweep_mode, fixed_angle
sweep_mode = filemetadata('sweep_mode')
fixed_angle = filemetadata('fixed_angle')
if scan_type == 'rhi':
sweep_mode['data'] = np.array(gfile.nsweeps * ['rhi'])
fixed_angle['data'] = gfile.how_attrs('azimuth', 'float32')
elif scan_type == 'ppi':
sweep_mode['data'] = np.array(gfile.nsweeps * ['azimuth_surveillance'])
fixed_angle['data'] = gfile.how_attrs('elevation', 'float32')
# time
time = filemetadata('time')
t_data = gfile.ray_header('timestamp', 'int64')
start_epoch = t_data[0] // 1.e6 # truncate to second resolution
start_time = datetime.datetime.utcfromtimestamp(start_epoch)
time['units'] = make_time_unit_str(start_time)
time['data'] = ((t_data - start_epoch * 1.e6) / 1.e6).astype('float64')
# range
_range = filemetadata('range')
ngates = int(gfile.raw_scan0_group_attr('how', 'bin_count'))
range_start = float(gfile.raw_scan0_group_attr('how', 'range_start'))
#n_samples insertion
range_samples = int(gfile.raw_scan0_group_attr('how', 'range_samples'))
range_step = float(gfile.raw_scan0_group_attr('how', 'range_step'))*range_samples
# range_step may need to be scaled by range_samples
# XXX This gives distances to start of gates not center, this matches
# Radx but may be incorrect, add range_step / 2. for center
_range['data'] = (np.arange(ngates, dtype='float32') * range_step +
range_start)
_range['meters_to_center_of_first_gate'] = range_start
_range['meters_between_gates'] = range_step
# elevation
elevation = filemetadata('elevation')
start_angle = gfile.ray_header('elevation_start', 'float32')
stop_angle = gfile.ray_header('elevation_stop', 'float32')
elevation['data'] = _avg_radial_angles(start_angle, stop_angle)
# azimuth
azimuth = filemetadata('azimuth')
start_angle = gfile.ray_header('azimuth_start', 'float32')
stop_angle = gfile.ray_header('azimuth_stop', 'float32')
azimuth['data'] = _avg_radial_angles(start_angle, stop_angle) % 360.
# fields
fields = {}
moment_groups = gfile.moment_groups()
moment_names = gfile.moment_names(moment_groups)
for moment_name, group in zip(moment_names, moment_groups):
field_name = filemetadata.get_field_name(moment_name)
if field_name is None:
continue
field_dic = filemetadata(field_name)
field_dic['data'] = gfile.moment_data(group, 'float32')
field_dic['_FillValue'] = get_fillvalue()
if valid_range_from_file:
try:
valid_min = gfile.raw_scan0_group_attr(group, 'dyn_range_min')
valid_max = gfile.raw_scan0_group_attr(group, 'dyn_range_max')
field_dic['valid_min'] = valid_min.decode('utf-8')
field_dic['valid_max'] = valid_max.decode('utf-8')
except:
pass
if units_from_file:
try:
units = gfile.raw_scan0_group_attr(group, 'unit')
field_dic['units'] = units.decode('utf-8')
except:
pass
fields[field_name] = field_dic
# ray_angle_res
ray_angle_res = filemetadata('ray_angle_res')
ray_angle_res['data'] = gfile.how_attrs('angle_step', 'float32')
# rays_are_indexed
rays_are_indexed = filemetadata('rays_are_indexed')
rays_are_indexed['data'] = np.array(
[['false', 'true'][i] for i in gfile.how_attrs('angle_sync', 'uint8')])
# target_scan_rate
target_scan_rate = filemetadata('target_scan_rate')
target_scan_rate['data'] = gfile.how_attrs('scan_speed', 'float32')
# scan_rate
scan_rate = filemetadata('scan_rate')
scan_rate['data'] = target_scan_rate['data']
if scan_type == 'ppi':
azs_names = ['az_speed', 'azimuth_speed']
azs_name = azs_names[0]
for azs_name in azs_names:
if gfile.is_field_in_ray_header(azs_name):
scan_rate['data'] = gfile.ray_header(azs_name, 'float32')
break
elif scan_type == 'rhi':
els_names = ['el_speed', 'elevation_speed']
els_name = els_names[0]
for els_name in els_names:
if gfile.is_field_in_ray_header(els_name):
scan_rate['data'] = gfile.ray_header(els_name, 'float32')
break
else:
scan_rate = None
# instrument_parameters
instrument_parameters = _get_instrument_params(gfile, filemetadata,
pulse_width)
gfile.close()
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,
ray_angle_res=ray_angle_res,
rays_are_indexed=rays_are_indexed,
scan_rate=scan_rate,
target_scan_rate=target_scan_rate)
def _get_instrument_params(gfile, filemetadata, pulse_width):
""" Return a dictionary containing instrument parameters. """
instrument_params = {}
dic = filemetadata('frequency')
dic['data'] = np.array(
[LIGHT_SPEED / gfile.raw_scan0_group_attr('how', 'radar_wave_length')],
dtype='float32')
instrument_params['frequency'] = dic
dic = filemetadata('radar_beam_width_h')
dic['data'] = gfile.how_attr('azimuth_beam', 'float32')
instrument_params['radar_beam_width_h'] = dic
dic = filemetadata('radar_beam_width_v')
dic['data'] = gfile.how_attr('elevation_beam', 'float32')
instrument_params['radar_beam_width_v'] = dic
dic = filemetadata('pulse_width')
pw_names = ['pulse_width_us', 'pulse_width_mks', 'pulse_width']
pw_name = 'pulse_width_us'
dic['data'] = None
for pw_name in pw_names:
if gfile.is_attr_in_group('/scan0/how', pw_name):
if pw_name == 'pulse_width':
dic['data'] = gfile.sweep_expand(
pulse_width[gfile.how_attrs(pw_name, 'int')[0]] * 1e-6)
else:
dic['data'] = gfile.sweep_expand(
gfile.how_attrs(pw_name, 'float32') * 1e-6)
break
instrument_params['pulse_width'] = dic
dic = filemetadata('prt')
dic['data'] = gfile.sweep_expand(1. / gfile.how_attrs('PRF', 'float32'))
instrument_params['prt'] = dic
unfolding = gfile.how_attrs('unfolding', 'int32')
dic = filemetadata('prt_mode')
dic['data'] = np.array([_prt_mode_from_unfolding(i) for i in unfolding])
instrument_params['prt_mode'] = dic
dic = filemetadata('prt_ratio')
dic['data'] = gfile.sweep_expand(
[[1, 2./3., 3./4., 4./5.][i] for i in unfolding])
instrument_params['prt_ratio'] = dic
dic = filemetadata('unambiguous_range')
dic['data'] = gfile.sweep_expand(gfile.how_attrs('range', 'float32'))
instrument_params['unambiguous_range'] = dic
if gfile.is_attr_in_group('/scan0/how/extended', 'nyquist_velocity'):
dic = filemetadata('nyquist_velocity')
dic['data'] = gfile.sweep_expand(
gfile.how_ext_attrs('nyquist_velocity'))
instrument_params['nyquist_velocity'] = dic
dic = filemetadata('n_samples')
dic['data'] = gfile.sweep_expand(
gfile.how_attrs('range_samples', 'int32') *
gfile.how_attrs('time_samples', 'int32'), dtype='int32')
instrument_params['n_samples'] = dic
return instrument_params
def _avg_radial_angles(angle1, angle2):
""" Return the average angle between two radial angles. """
return np.angle(
(np.exp(1.j*np.deg2rad(angle1)) +
np.exp(1.j*np.deg2rad(angle2))) / 2., deg=True)
def _prt_mode_from_unfolding(unfolding):
""" Return 'fixed' or 'staggered' depending on unfolding flag """
if unfolding == 0:
return 'fixed'
else:
return 'staggered'