<jupyter_start><jupyter_text>Test of vulipani phase processing<jupyter_code>import pyart
import numpy as np
from scipy import ndimage, interpolate
from matplotlib import pyplot as plt
%matplotlib inline
from netCDF4 import Dataset
input_file = '/home/rjackson/data/radar/cfrad.20110202_071003.000_to_20110202_071849.000_CPOL_PPI_level1b.nc'
radar = pyart.io.read(input_file)
radar.info()<jupyter_output>/home/rjackson/anaconda3/envs/lp_gatefiltertest/lib/python3.6/site-packages/jupyter_client/jsonutil.py:67: DeprecationWarning: Interpreting naive datetime as local 2018-08-24 15:02:58.688057. Please add timezone info to timestamps.
  new_obj[k] = extract_dates(v)
/home/rjackson/.local/lib/python3.6/site-packages/pyart/io/cfradial.py:376: RuntimeWarning: invalid value encountered in less
  data = self.ncvar[:]
/home/rjackson/.local/lib/python3.6/site-packages/pyart/io/cfradial.py:376: RuntimeWarning: invalid value encountered in greater
  data = self.ncvar[:]
/home/rjackson/.local/lib/python3.6/site-packages/pyart/io/cfradial.py:376: UserWarning: WARNING: valid_min not used since it
cannot be safely cast to variable data type
  data = self.ncvar[:]
/home/rjackson/.local/lib/python3.6/site-packages/pyart/io/cfradial.py:376: UserWarning: WARNING: valid_max not used since it
cannot be safely cast to variable data type
  data = self.ncvar[:]<jupyter_text>Derive pseudo-NCP<jupyter_code>print('Calculating texture....')
nyq_Gunn = radar.instrument_parameters['nyquist_velocity']['data'][0]
data = ndimage.filters.generic_filter(radar.fields['velocity']['data'],
                                      pyart.util.interval_std, size = (3,3),
                                      extra_arguments = (-nyq_Gunn, nyq_Gunn))
filtered_data = ndimage.filters.median_filter(data, size = (3,3))
texture_field = pyart.config.get_metadata('corrected_velocity')
texture_field['data'] = filtered_data
texture_field['data'] = filtered_data
radar.add_field('velocity_texture', texture_field, replace_existing = True)
texture = radar.fields['velocity_texture']['data']
plt.hist(texture[texture > 0], normed=1, bins=40)
plt.xlabel('Velocity texture')
plt.ylabel('p.d.f')
## For testing, assume zero psi_dp field
psidp = np.zeros(texture.shape)
radar.add_field_like('PHIDP', 'PSIDP', NCP, replace_existing = True)
NCP = np.zeros(texture.shape)
NCP[texture < 2] = 1
radar.add_field_like('velocity_texture', 'normalized_coherent_power', NCP, replace_existing = True)
display = pyart.graph.RadarDisplay(radar)
plt.figure(figsize=(6,6))
display.plot_ppi('bringi_differential_phase', vmin=0, vmax=360, cmap=pyart.graph.cm.NWSRef, gatefilter=gatefilter)
plt.xlim([-50,50])
plt.ylim([-50,50])<jupyter_output>/home/rjackson/anaconda3/envs/lp_gatefiltertest/lib/python3.6/site-packages/jupyter_client/jsonutil.py:67: DeprecationWarning: Interpreting naive datetime as local 2018-08-24 15:01:36.746285. Please add timezone info to timestamps.
  new_obj[k] = extract_dates(v)<jupyter_text>Derive KDP<jupyter_code>gatefilter = pyart.filters.GateFilter(radar)
gatefilter.exclude_above('velocity_texture', 2)
gatefilter.exclude_below('cross_correlation_ratio', 0.9)
gatefilter.exclude_masked('corrected_reflectivity')
gatefilter = pyart.correct.despeckle_field(radar, 'corrected_reflectivity', gatefilter=gatefilter, size=30)
#filter_dict = {}
#filter_dict['rhohv_field'] = 'RHOHV'


kdp_dict = pyart.correct.phase_proc_lp_gf(radar, gatefilter=gatefilter, min_rhv=0.95,
                              LP_solver='cylp', refl_field='corrected_reflectivity')
print(kdp_dict)
radar.add_field('specific_differential_phase', kdp_dict[1], replace_existing = True)
radar.add_field('corrected_differential_phase', kdp_dict[0], replace_existing = True)
plt.figure(figsize=(6,13))
plt.subplot(212)
display = pyart.graph.RadarDisplay(radar)
display.plot_ppi('corrected_differential_phase', vmin=-90, vmax=90, gatefilter=gatefilter,)
plt.xlim([-50,50])
plt.ylim([-50,50])

plt.subplot(211)
display.plot_ppi('specific_differential_phase', vmin=0, vmax=1, gatefilter=gatefilter,
                 cmap=pyart.graph.cm_colorblind.HomeyerRainbow)
plt.xlim([-50,50])
plt.ylim([-50,50])
altitudes = radar.gate_altitude['data']

sounding_input_file = '/home/rjackson/data/soundings/twpsondewnpnC3.b1.20111202.111900.custom.cdf'
Sounding_netcdf = Dataset(sounding_input_file, mode='r')

# Convert timestamps to datetime format
Time = Sounding_netcdf.variables['time_offset'][:]
base_time = Sounding_netcdf.variables['base_time'][:]
alt = Sounding_netcdf.variables['alt'][:]
u = Sounding_netcdf.variables['u_wind'][:]
v = Sounding_netcdf.variables['v_wind'][:]
t = Sounding_netcdf.variables['tdry'][:]

function = interpolate.interp1d(alt,t,bounds_error=False)
temp = np.zeros(altitudes.shape)

for j in range(0,592):
    temp[:,j] = function(altitudes[:,j])

radar.add_field_like('specific_differential_phase', 
                     'temperature', temp, replace_existing = True)
display.plot_ppi('temperature', vmin=-20, vmax=30, sweep=2,
                 cmap=pyart.graph.cm.NWSRef)
#plt.plot(alt,t)
hydrometeor_class = pyart.retrieve.hydroclass_semisupervised(radar, 
                                                             refl_field='Refl',
                                                             rhv_field='RHOHV',
                                                             temp_field='temperature',
                                                             zdr_field='ZDR',
                                                             kdp_field='corrected_differential_phase',
                                                             weights=np.array([1.0, 
                                                                               0.75,
                                                                               1.,
                                                                               0.75,
                                                                               0.5]))
radar.add_field('hydrometeor_class', hydrometeor_class, replace_existing=True)
radar_xsection = pyart.util.cross_section_ppi(radar, [200])

display_xsection = pyart.graph.RadarMapDisplay(radar_xsection)
plt.figure(figsize=(6,12))
plt.subplot(211)
display_xsection.plot_rhi('hydrometeor_class',
                          vmin=0,
                          vmax=8)
plt.ylim([0,20])
plt.xlim([0,80])
print(np.max(hydrometeor_class))
display_xsection.plot_colorbar(ticks=[0,1,2,3,4,5,6,7,8],
                               ticklabs=['DS','CR', 
                                         'LR','GR',
                                         'RN','VI',
                                         'WS','MH',
                                         'IH'])

plt.subplot(212)
display_xsection.plot_rhi('Refl')
plt.ylim([0,20])
plt.xlim([0,80])
print(hydrometeor_class)<jupyter_output>{'units': 'legend', 'long_name': 'Radar Echo classification', 'standard_name': 'radar_echo_classification', 'data': array([[0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       ..., 
       [3, 3, 0, ..., 0, 0, 0],
       [3, 3, 0, ..., 0, 0, 0],
       [3, 3, 0, ..., 0, 0, 0]]), 'coordinates': 'elevation azimuth range'}
{'units': 'legend', 'long_name': 'Radar Echo classification', 'standard_name': 'radar_echo_classification', 'data': array([[0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       [0, 0, 0, ..., 0, 0, 0],
       ..., 
       [3, 3, 0, ..., 0, 0, 0],
       [3, 3, 0, ..., 0, 0, 0],
       [3, 3, 0, ..., 0, 0, 0]]), 'coordinates': 'elevation azimuth range'}
{'fraction_non_alphanumeric': 0.14048059149722736, 'fraction_numerical': 0.04172167942962767, 'mean_word_length': 2.868743615934627, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 11, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 34, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '3', 'dir': 'jupyter-structured-clean-dedup', 'id': '394311', 'license': 'MIT', 'n_tokens_mistral': 2803, 'n_tokens_neox': 2503, 'n_words': 417, 'path': 'notebooks/Phase processing vulipani.ipynb', 'repo_name': 'EVS-ATMOS/cmdv-rrm-anl', 'size': '7574'}