diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index 86b22e1ffa..fcb51b3349 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -104,6 +104,7 @@ def to_vpt(radar, single_scan=True): return + def join_radar(radar1, radar2): """ Combine two radar instances into one. @@ -135,6 +136,38 @@ def join_radar(radar1, radar2): new_radar.nsweeps += radar2.nsweeps new_radar.sweep_mode['data'] = np.append( radar1.sweep_mode['data'], radar2.sweep_mode['data']) + if ((radar1.rays_are_indexed is not None) and + (radar2.rays_are_indexed is not None)): + new_radar.rays_are_indexed['data'] = np.append( + radar1.rays_are_indexed['data'], + radar2.rays_are_indexed['data']) + else: + new_radar.rays_are_indexed = None + + if new_radar.instrument_parameters is not None: + if 'nyquist_velocity' in new_radar.instrument_parameters: + new_radar.instrument_parameters['nyquist_velocity']['data'] = ( + np.append( + radar1.instrument_parameters['nyquist_velocity']['data'], + radar2.instrument_parameters['nyquist_velocity']['data'])) + if 'pulse_width' in new_radar.instrument_parameters: + new_radar.instrument_parameters['pulse_width']['data'] = ( + np.append( + radar1.instrument_parameters['pulse_width']['data'], + radar2.instrument_parameters['pulse_width']['data'])) + if 'number_of_pulses' in new_radar.instrument_parameters: + new_radar.instrument_parameters['number_of_pulses']['data'] = ( + np.append( + radar1.instrument_parameters['number_of_pulses']['data'], + radar2.instrument_parameters['number_of_pulses']['data'])) + + if ((radar1.ray_angle_res is not None) and + (radar2.ray_angle_res is not None)): + new_radar.ray_angle_res['data'] = np.append( + radar1.ray_angle_res['data'], + radar2.ray_angle_res['data']) + else: + new_radar.ray_angle_res = None if len(radar1.range['data']) >= len(radar2.range['data']): new_radar.range['data'] = radar1.range['data'] @@ -151,16 +184,26 @@ def join_radar(radar1, radar2): new_radar.time['units'] = datetime_utils.EPOCH_UNITS new_radar.nrays = len(new_radar.time['data']) + fields_to_remove = [] for var in new_radar.fields.keys(): - sh1 = radar1.fields[var]['data'].shape - sh2 = radar2.fields[var]['data'].shape - new_field_shape = (sh1[0] + sh2[0], max(sh1[1], sh2[1])) - new_field = np.ma.zeros(new_field_shape) - new_field[:] = np.ma.masked - new_field.set_fill_value(get_fillvalue()) - new_field[0:sh1[0], 0:sh1[1]] = radar1.fields[var]['data'] - new_field[sh1[0]:, 0:sh2[1]] = radar2.fields[var]['data'] - new_radar.fields[var]['data'] = new_field + # if the field is present in both radars combine both fields + # otherwise remove it from new radar + if var in radar1.fields and var in radar2.fields: + sh1 = radar1.fields[var]['data'].shape + sh2 = radar2.fields[var]['data'].shape + new_field_shape = (sh1[0] + sh2[0], max(sh1[1], sh2[1])) + new_field = np.ma.masked_all(new_field_shape) + new_field.set_fill_value(get_fillvalue()) + new_field[0:sh1[0], 0:sh1[1]] = radar1.fields[var]['data'] + new_field[sh1[0]:, 0:sh2[1]] = radar2.fields[var]['data'] + new_radar.fields[var]['data'] = new_field + else: + print("Field "+var+" not present in both radars") + fields_to_remove.append(var) + + if fields_to_remove: + for field_name in fields_to_remove: + new_radar.fields.pop(field_name, None) # radar locations # TODO moving platforms - any more? @@ -196,8 +239,10 @@ def join_radar(radar1, radar2): radar2.longitude['data']) new_radar.altitude['data'] = np.append(radar1.altitude['data'], radar2.altitude['data']) + return new_radar + def image_mute_radar( radar, field, mute_field, mute_threshold, field_threshold=None): """ diff --git a/pyart/util/tests/test_radar_utils.py b/pyart/util/tests/test_radar_utils.py index 5a38537cf0..e3bf52bc2f 100644 --- a/pyart/util/tests/test_radar_utils.py +++ b/pyart/util/tests/test_radar_utils.py @@ -10,6 +10,24 @@ def test_is_vpt(): pyart.util.to_vpt(radar) assert pyart.util.is_vpt(radar) +def test_join_radar(): + radar1 = pyart.testing.make_empty_ppi_radar(10, 36, 3) + radar1.instrument_parameters = { + 'nyquist_velocity': {'data': np.array([6] * 3)} + } + field = {'data': np.ones((36 * 3, 10))} + radar1.add_field('f1', field) + radar1.add_field('f2', field) + radar2 = pyart.testing.make_empty_ppi_radar(10, 36, 4) + radar2.instrument_parameters = { + 'nyquist_velocity': {'data': np.array([8] * 4)} + } + field = {'data': np.ones((36 * 4, 10))} + radar2.add_field('f1', field) + + radar3 = pyart.util.join_radar(radar1, radar2) + assert 'f1' in radar3.fields + assert len(radar3.instrument_parameters['nyquist_velocity']['data']) == 7 def test_to_vpt(): # single scan