"""
Utilities for working with phase data.
Code based upon algorithm descriped in:
S. E. Giangrande et al, J. of Atmos. and Ocean. Tech., 2013, 30, 1716.
Adapted by Scott Collis and Scott Giangrande, refactored by Jonathan Helmus.
"""
import copy
from time import time
import numpy as np
import scipy.ndimage
from numpy import ma
from ..config import get_field_name, get_fillvalue, get_metadata
from ..filters import GateFilter
from ..util import rolling_window
def det_sys_phase(
radar, ncp_lev=0.4, rhohv_lev=0.6, ncp_field=None, rhv_field=None, phidp_field=None
):
"""
Determine the system phase.
Parameters
----------
radar : Radar
Radar object for which to determine the system phase.
ncp_lev : float, optional
Miminum normal coherent power level. Regions below this value will
not be included in the phase calculation.
rhohv_lev : float, optional
Miminum copolar coefficient level. Regions below this value will not
be included in the phase calculation.
ncp_field, rhv_field, phidp_field : str, optional
Field names within the radar object which represent the normal
coherent power, the copolar coefficient, and the differential phase
shift. A value of None for any of these parameters will use the
default field name as defined in the Py-ART configuration file.
Returns
-------
sys_phase : float or None
Estimate of the system phase. None is not estimate can be made.
"""
# parse the field parameters
if ncp_field is None:
ncp_field = get_field_name("normalized_coherent_power")
if rhv_field is None:
rhv_field = get_field_name("cross_correlation_ratio")
if phidp_field is None:
phidp_field = get_field_name("differential_phase")
ncp = radar.fields[ncp_field]["data"][:, 30:]
rhv = radar.fields[rhv_field]["data"][:, 30:]
phidp = radar.fields[phidp_field]["data"][:, 30:]
last_ray_idx = radar.sweep_end_ray_index["data"][0]
return _det_sys_phase(ncp, rhv, phidp, last_ray_idx, ncp_lev, rhohv_lev)
def _det_sys_phase(ncp, rhv, phidp, last_ray_idx, ncp_lev=0.4, rhv_lev=0.6):
"""Determine the system phase, see :py:func:`det_sys_phase`."""
good = False
phases = []
for radial in range(last_ray_idx + 1):
meteo = np.logical_and(ncp[radial, :] > ncp_lev, rhv[radial, :] > rhv_lev)
mpts = np.where(meteo)
if len(mpts[0]) > 25:
good = True
msmth_phidp = smooth_and_trim(phidp[radial, mpts[0]], 9)
phases.append(msmth_phidp[0:25].min())
if not good:
return None
return np.median(phases)
def fzl_index(fzl, ranges, elevation, radar_height):
"""
Return the index of the last gate below a given altitude.
Parameters
----------
fzl : float
Maximum altitude.
ranges : array
Range to measurement volume/gate in meters.
elevation : float
Elevation of antenna in degrees.
radar_height :
Altitude of radar in meters.
Returns
-------
idx : int
Index of last gate which has an altitude below `fzl`.
Notes
-----
Standard atmosphere is assumed, R = 4 / 3 * Re
"""
Re = 6371.0 * 1000.0
p_r = 4.0 * Re / 3.0
z = (
radar_height
+ (
ranges**2
+ p_r**2
+ 2.0 * ranges * p_r * np.sin(elevation * np.pi / 180.0)
)
** 0.5
- p_r
)
# Make sure the freezing level isn't under the radar!
# Return the minimum window size for the 5-pt filter
if np.all(z > fzl):
return 6
else:
return np.where(z < fzl)[0].max()
def det_process_range(radar, sweep, fzl, doc=10):
"""
Determine the processing range for a given sweep.
Queues the radar and returns the indices which can be used to slice
the radar fields and select the desired sweep with gates which are
below a given altitude.
Parameters
----------
radar : Radar
Radar object from which ranges will be determined.
sweep : int
Sweep (0 indexed) for which to determine processing ranges.
fzl : float
Maximum altitude in meters. The determined range will not include
gates which are above this limit.
doc : int, optional
Minimum number of gates which will be excluded from the determined
range.
Returns
-------
gate_end : int
Index of last gate below `fzl` and satisfying the `doc` parameter.
ray_start : int
Ray index which defines the start of the region.
ray_end : int
Ray index which defined the end of the region.
"""
# determine the index of the last valid gate
ranges = radar.range["data"]
elevation = radar.fixed_angle["data"][sweep]
radar_height = radar.altitude["data"]
gate_end = fzl_index(fzl, ranges, elevation, radar_height)
if doc is not None:
gate_end = min(gate_end, len(ranges) - doc)
else:
gate_end = min(gate_end, len(ranges))
ray_start = radar.sweep_start_ray_index["data"][sweep]
ray_end = radar.sweep_end_ray_index["data"][sweep] + 1
return gate_end, ray_start, ray_end
def snr(line, wl=11):
"""Return the signal to noise ratio after smoothing."""
signal = smooth_and_trim(line, window_len=wl)
_noise = smooth_and_trim(np.sqrt((line - signal) ** 2), window_len=wl)
return abs(signal) / _noise
def smooth_masked(raw_data, wind_len=11, min_valid=6, wind_type="median"):
"""
Smoothes the data using a rolling window.
data with less than n valid points is masked.
Parameters
----------
raw_data : float masked array
The data to smooth.
win_len : float
Length of the moving window.
min_valid : float
Minimum number of valid points for the smoothing to be valid.
wind_type : str
Type of window. Can be median or mean.
Returns
-------
data_smooth : float masked array
Smoothed data.
"""
valid_wind = ["median", "mean"]
if wind_type not in valid_wind:
raise ValueError("Window " + wind_type + " is none of " + " ".join(valid_wind))
# we want an odd window
if wind_len % 2 == 0:
wind_len += 1
half_wind = int((wind_len - 1) / 2)
# initialize smoothed data
nrays, nbins = np.shape(raw_data)
data_smooth = np.ma.zeros((nrays, nbins))
data_smooth[:] = np.ma.masked
data_smooth.set_fill_value(get_fillvalue())
mask = np.ma.getmaskarray(raw_data)
valid = np.logical_not(mask)
mask_wind = rolling_window(mask, wind_len)
valid_wind = np.logical_not(mask_wind).astype(int)
nvalid = np.sum(valid_wind, -1)
data_wind = rolling_window(raw_data, wind_len)
# check which gates are valid
ind_valid = np.logical_and(
nvalid >= min_valid, valid[:, half_wind:-half_wind]
).nonzero()
if data_wind is not None:
data_smooth[ind_valid[0], ind_valid[1] + half_wind] = eval(
"np.ma." + wind_type + "(data_wind, axis=-1)"
)[ind_valid]
return data_smooth
def unwrap_masked(lon, centered=False, copy=True):
"""
Unwrap a sequence of longitudes or headings in degrees.
Parameters
----------
lon : array
Longtiudes or heading in degress. If masked output will also be
masked.
centered : bool, optional
Center the unwrapping as close to zero as possible.
copy : bool, optional.
True to return a copy, False will avoid a copy when possible.
Returns
-------
unwrap : array
Array of unwrapped longtitudes or headings, in degrees.
"""
masked_input = ma.isMaskedArray(lon)
if masked_input:
fill_value = lon.fill_value
# masked_invalid loses the original fill_value (ma bug, 2011/01/20)
lon = np.ma.masked_invalid(lon).astype(float)
if lon.ndim != 1:
raise ValueError("Only 1-D sequences are supported")
if lon.shape[0] < 2:
return lon
x = lon.compressed()
if len(x) < 2:
return lon
w = np.zeros(x.shape[0] - 1, int)
ld = np.diff(x)
np.putmask(w, ld > 180, -1)
np.putmask(w, ld < -180, 1)
x[1:] += w.cumsum() * 360.0
if centered:
x -= 360 * np.round(x.mean() / 360.0)
if lon.mask is ma.nomask:
lon[:] = x
else:
lon[~lon.mask] = x
if masked_input:
lon.fill_value = fill_value
return lon
else:
return lon.filled(np.nan)
# this function adapted from the Scipy Cookbook:
# http://www.scipy.org/Cookbook/SignalSmooth
def smooth_and_trim(x, window_len=11, window="hanning"):
"""
Smooth data using a window with requested size.
This method is based on the convolution of a scaled window with the signal.
The signal is prepared by introducing reflected copies of the signal
(with the window size) in both ends so that transient parts are minimized
in the begining and end part of the output signal.
Parameters
----------
x : array
The input signal.
window_len : int, optional
The dimension of the smoothing window; should be an odd integer.
window : str
The type of window from 'flat', 'hanning', 'hamming', 'bartlett',
'blackman' or 'sg_smooth'. A flat window will produce a moving
average smoothing.
Returns
-------
y : array
The smoothed signal with length equal to the input signal.
"""
if x.ndim != 1:
raise ValueError("smooth only accepts 1 dimension arrays.")
if x.size < window_len:
raise ValueError("Input vector needs to be bigger than window size.")
if window_len < 3:
return x
valid_windows = ["flat", "hanning", "hamming", "bartlett", "blackman", "sg_smooth"]
if window not in valid_windows:
raise ValueError("Window is on of " + " ".join(valid_windows))
s = np.r_[x[window_len - 1 : 0 : -1], x, x[-1:-window_len:-1]]
if window == "flat": # moving average
w = np.ones(int(window_len), "d")
elif window == "sg_smooth":
w = np.array([0.1, 0.25, 0.3, 0.25, 0.1])
else:
w = eval("np." + window + "(window_len)")
y = np.convolve(w / w.sum(), s, mode="valid")
return y[int(window_len / 2) : len(x) + int(window_len / 2)]
# adapted smooth and trim function to work with 2dimensional arrays
def smooth_and_trim_scan(x, window_len=11, window="hanning"):
"""
Smooth data using a window with requested size.
This method is based on the convolution of a scaled window with the signal.
The signal is prepared by introducing reflected copies of the signal
(with the window size) in both ends so that transient parts are minimized
in the begining and end part of the output signal.
Parameters
----------
x : ndarray
The input signal.
window_len : int, optional
The dimension of the smoothing window; should be an odd integer.
window : str, optional
The type of window from 'flat', 'hanning', 'hamming', 'bartlett',
'blackman' or 'sg_smooth'. A flat window will produce a moving
average smoothing.
Returns
-------
y : ndarray
The smoothed signal with length equal to the input signal.
"""
from scipy.ndimage import convolve1d
if x.ndim != 2:
raise ValueError("smooth only accepts 2 dimension arrays.")
if x.shape[1] < window_len:
mess = "Input dimension 1 needs to be bigger than window size."
raise ValueError(mess)
if window_len < 3:
return x
valid_windows = ["flat", "hanning", "hamming", "bartlett", "blackman", "sg_smooth"]
if window not in valid_windows:
raise ValueError("Window is on of " + " ".join(valid_windows))
if window == "flat": # moving average
w = np.ones(int(window_len), "d")
elif window == "sg_smooth":
w = np.array([0.1, 0.25, 0.3, 0.25, 0.1])
else:
w = eval("np." + window + "(window_len)")
y = convolve1d(x, w / w.sum(), axis=1)
return y
def noise(line, wl=11):
"""Return the noise after smoothing."""
signal = smooth_and_trim(line, window_len=wl)
_noise = np.sqrt((line - signal) ** 2)
return _noise
def get_phidp_unf(
radar,
ncp_lev=0.4,
rhohv_lev=0.6,
debug=False,
ncpts=20,
doc=-10,
overide_sys_phase=False,
sys_phase=-135,
nowrap=None,
refl_field=None,
ncp_field=None,
rhv_field=None,
phidp_field=None,
):
"""
Get Unfolded Phi differential phase
Parameters
----------
radar : Radar
The input radar.
ncp_lev : float, optional
Miminum normal coherent power level. Regions below this value will
not be included in the calculation.
rhohv_lev : float, optional
Miminum copolar coefficient level. Regions below this value will not
be included in the calculation.
debug : bool, optional
True to print debugging information, False to supress printing.
ncpts : int, optional
Minimum number of points in a ray. Regions within a ray smaller than
this or beginning before this gate number are excluded from
calculations.
doc : int or None, optional
Index of first gate not to include in field data, None include all.
overide_sys_phase : bool, optional
True to use `sys_phase` as the system phase. False will determine a
value automatically.
sys_phase : float, optional
System phase, not used if overide_sys_phase is False.
nowrap : int or None, optional
Gate number where unwrapping should begin. `None` will unwrap all
gates.
refl_field ncp_field, rhv_field, phidp_field : str, optional
Field names within the radar object which represent the horizonal
reflectivity, normal coherent power, the copolar coefficient, and the
differential phase shift. A value of None for any of these parameters
will use the default field name as defined in the Py-ART
configuration file.
Returns
-------
cordata : array
Unwrapped phi differential phase.
"""
# parse the field parameters
if refl_field is None:
refl_field = get_field_name("reflectivity")
if ncp_field is None:
ncp_field = get_field_name("normalized_coherent_power")
if rhv_field is None:
rhv_field = get_field_name("cross_correlation_ratio")
if phidp_field is None:
phidp_field = get_field_name("differential_phase")
if doc is not None:
my_phidp = radar.fields[phidp_field]["data"][:, 0:doc]
my_rhv = radar.fields[rhv_field]["data"][:, 0:doc]
my_ncp = radar.fields[ncp_field]["data"][:, 0:doc]
my_z = radar.fields[refl_field]["data"][:, 0:doc]
else:
my_phidp = radar.fields[phidp_field]["data"]
my_rhv = radar.fields[rhv_field]["data"]
my_ncp = radar.fields[ncp_field]["data"]
my_z = radar.fields[refl_field]["data"]
t = time()
if overide_sys_phase:
system_zero = sys_phase
else:
system_zero = det_sys_phase(
radar, ncp_field=ncp_field, rhv_field=rhv_field, phidp_field=phidp_field
)
if system_zero is None:
system_zero = sys_phase
cordata = np.zeros(my_rhv.shape, dtype=float)
for radial in range(my_rhv.shape[0]):
my_snr = snr(my_z[radial, :])
notmeteo = np.logical_or(
np.logical_or(my_ncp[radial, :] < ncp_lev, my_rhv[radial, :] < rhohv_lev),
my_snr < 10.0,
)
x_ma = ma.masked_where(notmeteo, my_phidp[radial, :])
try:
ma.notmasked_contiguous(x_ma)
for slc in ma.notmasked_contiguous(x_ma):
# so trying to get rid of clutter and small things that
# should not add to phidp anyway
if slc.stop - slc.start < ncpts or slc.start < ncpts:
x_ma.mask[slc.start - 1 : slc.stop + 1] = True
c = 0
except TypeError: # non sequence, no valid regions
c = 1 # ie do nothing
x_ma.mask = True
except AttributeError:
# sys.stderr.write('No Valid Regions, ATTERR \n ')
# sys.stderr.write(myfile.times['time_end'].strftime('%Y-%m-%dT%H:%M:%SZ') + '\n')
# print x_ma
# print x_ma.mask
c = 1 # also do nothing
x_ma.mask = True
if nowrap is not None:
# Start the unfolding a bit later in order to avoid false
# jumps based on clutter
unwrapped = copy.deepcopy(x_ma)
end_unwrap = unwrap_masked(x_ma[nowrap::], centered=False)
unwrapped[nowrap::] = end_unwrap
else:
unwrapped = unwrap_masked(x_ma, centered=False)
# end so no clutter expected
system_max = (
unwrapped[np.where(np.logical_not(notmeteo))][-10:-1].mean() - system_zero
)
unwrapped_fixed = np.zeros(len(x_ma), dtype=float)
based = unwrapped - system_zero
based[0] = 0.0
notmeteo[0] = False
based[-1] = system_max
notmeteo[-1] = False
unwrapped_fixed[np.where(np.logical_not(based.mask))[0]] = based[
np.where(np.logical_not(based.mask))[0]
]
if len(based[np.where(np.logical_not(based.mask))[0]]) > 11:
unwrapped_fixed[np.where(based.mask)[0]] = np.interp(
np.where(based.mask)[0],
np.where(np.logical_not(based.mask))[0],
smooth_and_trim(based[np.where(np.logical_not(based.mask))[0]]),
)
else:
unwrapped_fixed[np.where(based.mask)[0]] = np.interp(
np.where(based.mask)[0],
np.where(np.logical_not(based.mask))[0],
based[np.where(np.logical_not(based.mask))[0]],
)
if c != 1:
cordata[radial, :] = unwrapped_fixed
else:
cordata[radial, :] = np.zeros(my_rhv.shape[1])
if debug:
print("Exec time: ", time() - t)
return cordata
def construct_A_matrix(n_gates, filt):
"""
Construct a row-augmented A matrix. Equation 5 in Giangrande et al, 2012.
A is a block matrix given by:
.. math::
\\bf{A} = \\begin{bmatrix} \\bf{I} & \\bf{-I} \\\\\\\\
\\bf{-I} & \\bf{I} \\\\\\\\ \\bf{Z}
& \\bf{M} \\end{bmatrix}
where
:math:`\\bf{I}` is the identity matrix
:math:`\\bf{Z}` is a matrix of zeros
:math:`\\bf{M}` contains our differential constraints.
Each block is of shape n_gates by n_gates making
shape(:math:`\\bf{A}`) = (3 * n, 2 * n).
Note that :math:`\\bf{M}` contains some side padding to deal with edge
issues.
Parameters
----------
n_gates : int
Number of gates, determines size of identity matrix.
filt : array
Input filter.
Returns
-------
a : matrix
Row-augmented A matrix.
"""
Identity = np.eye(n_gates)
filter_length = len(filt)
M_matrix_middle = np.diag(np.ones(n_gates - filter_length + 1), k=0) * 0.0
posn = np.linspace(
-1.0 * (filter_length - 1) / 2, (filter_length - 1) / 2, filter_length
)
for diag in range(filter_length):
M_matrix_middle = (
M_matrix_middle
+ np.diag(
np.ones(int(n_gates - filter_length + 1 - np.abs(posn[diag]))),
k=int(posn[diag]),
)
* filt[diag]
)
side_pad = (filter_length - 1) // 2
M_matrix = np.bmat(
[
np.zeros([n_gates - filter_length + 1, side_pad], dtype=float),
M_matrix_middle,
np.zeros([n_gates - filter_length + 1, side_pad], dtype=float),
]
)
Z_matrix = np.zeros([n_gates - filter_length + 1, n_gates])
return np.bmat(
[[Identity, -1.0 * Identity], [Identity, Identity], [Z_matrix, M_matrix]]
)
def construct_B_vectors(phidp_mod, z_mod, filt, coef=0.914, dweight=60000.0):
"""
Construct B vectors. See Giangrande et al, 2012.
Parameters
----------
phidp_mod : 2D array
Phi differential phases.
z_mod : 2D array.
Reflectivity, modified as needed.
filt : array
Input filter.
coef : float, optional.
Cost coefficients.
dweight : float, optional.
Weights.
Returns
-------
b : matrix
Matrix containing B vectors.
"""
n_gates = phidp_mod.shape[1]
n_rays = phidp_mod.shape[0]
filter_length = len(filt)
side_pad = (filter_length - 1) // 2
top_of_B_vectors = np.bmat([[-phidp_mod, phidp_mod]])
data_edges = np.bmat(
[
phidp_mod[:, 0:side_pad],
np.zeros([n_rays, n_gates - filter_length + 1]),
phidp_mod[:, -side_pad:],
]
)
ii = filter_length - 1
jj = data_edges.shape[1] - 1
list_corrl = np.zeros([n_rays, jj - ii + 1])
for count in range(list_corrl.shape[1]):
list_corrl[:, count] = -1.0 * (
np.array(filt) * (np.asarray(data_edges))[:, count : count + ii + 1]
).sum(axis=1)
sct = ((10.0 ** (0.1 * z_mod)) ** coef / dweight)[:, side_pad:-side_pad]
sct[np.where(sct < 0.0)] = 0.0
sct[:, 0:side_pad] = list_corrl[:, 0:side_pad]
sct[:, -side_pad:] = list_corrl[:, -side_pad:]
B_vectors = np.bmat([[top_of_B_vectors, sct]])
return B_vectors
def LP_solver_cvxopt(A_Matrix, B_vectors, weights, solver="glpk"):
"""
Solve the Linear Programming problem given in Giangrande et al, 2012 using
the CVXOPT module.
Parameters
----------
A_Matrix : matrix
Row augmented A matrix, see :py:func:`construct_A_matrix`
B_vectors : matrix
Matrix containing B vectors, see :py:func:`construct_B_vectors`
weights : array
Weights.
solver : str or None
LP solver backend to use, choices are 'glpk', 'mosek' or None to use
the conelp function in CVXOPT. 'glpk' and 'mosek' are only available
if they are installed and CVXOPT was build with the correct bindings.
Returns
-------
soln : array
Solution to LP problem.
See Also
--------
LP_solver_pyglpk : Solve LP problem using the PyGLPK module.
LP_solver_cylp : Solve LP problem using the cylp module.
LP_solver_cylp_mp : Solve LP problem using the cylp module
using multi processes.
"""
from cvxopt import matrix, solvers
n_gates = weights.shape[1] // 2
n_rays = B_vectors.shape[0]
mysoln = np.zeros([n_rays, n_gates])
G = matrix(np.bmat([[-A_Matrix], [-np.eye(2 * n_gates)]]))
h_array = np.zeros(5 * n_gates - 4)
for raynum in range(n_rays):
c = matrix(weights[raynum]).T
h_array[: 3 * n_gates - 4] = -B_vectors[raynum]
h = matrix(h_array)
sol = solvers.lp(c, G, h, solver=solver)
# XXX when a solution is not found sol is None, need to check and
# deal with this...
# extract the solution
this_soln = np.zeros(n_gates)
for i in range(n_gates):
this_soln[i] = sol["x"][i + n_gates]
# apply smoothing filter and record in output array
mysoln[raynum, :] = smooth_and_trim(this_soln, window_len=5, window="sg_smooth")
return mysoln
def LP_solver_pyglpk(
A_Matrix, B_vectors, weights, it_lim=7000, presolve=True, really_verbose=False
):
"""
Solve the Linear Programming problem given in Giangrande et al, 2012 using
the PyGLPK module.
Parameters
----------
A_Matrix : matrix
Row augmented A matrix, see :py:func:`construct_A_matrix`
B_vectors : matrix
Matrix containing B vectors, see :py:func:`construct_B_vectors`
weights : array
Weights.
it_lim : int, optional
Simplex iteration limit.
presolve : bool, optional
True to use the LP presolver.
really_verbose : bool, optional
True to print LPX messaging. False to suppress.
Returns
-------
soln : array
Solution to LP problem.
See Also
--------
LP_solver_cvxopt : Solve LP problem using the CVXOPT module.
LP_solver_cylp : Solve LP problem using the cylp module.
LP_solver_cylp_mp : Solve LP problem using the cylp module
using multi processes.
"""
import glpk
if really_verbose:
message_state = glpk.LPX.MSG_ON
else:
message_state = glpk.LPX.MSG_OFF
n_gates = weights.shape[1] // 2
n_rays = B_vectors.shape[0]
mysoln = np.zeros([n_rays, n_gates])
lp = glpk.LPX() # Create empty problem instance
lp.name = "LP_MIN" # Assign symbolic name to problem
lp.obj.maximize = False # Set this as a maximization problem
lp.rows.add(2 * n_gates + n_gates - 4) # Append rows
lp.cols.add(2 * n_gates)
glpk.env.term_on = True
for cur_row in range(2 * n_gates + n_gates - 4):
lp.rows[cur_row].matrix = list(np.squeeze(np.asarray(A_Matrix[cur_row, :])))
for i in range(2 * n_gates):
lp.cols[i].bounds = 0.0, None
for raynum in range(n_rays):
this_soln = np.zeros(n_gates)
for i in range(2 * n_gates + n_gates - 4):
lp.rows[i].bounds = B_vectors[raynum, i], None
for i in range(2 * n_gates):
lp.obj[i] = weights[raynum, i]
lp.simplex(
msg_lev=message_state,
meth=glpk.LPX.PRIMAL,
it_lim=it_lim,
presolve=presolve,
)
for i in range(n_gates):
this_soln[i] = lp.cols[i + n_gates].primal
mysoln[raynum, :] = smooth_and_trim(this_soln, window_len=5, window="sg_smooth")
return mysoln
def solve_cylp(model, B_vectors, weights, ray, chunksize):
"""
Worker process for LP_solver_cylp_mp.
Parameters
----------
model : CyClpModel
Model of the LP Problem, see :py:func:`LP_solver_cylp_mp`
B_vectors : matrix
Matrix containing B vectors, see :py:func:`construct_B_vectors`
weights : array
Weights.
ray : int
Starting ray.
chunksize : int
Number of rays to process.
Returns
-------
soln : array
Solution to LP problem.
See Also
--------
LP_solver_cylp_mp : Parent function.
LP_solver_cylp : Single Process Solver.
"""
from cylp.cy.CyClpSimplex import CyClpSimplex
n_gates = weights.shape[1] // 2
soln = np.zeros([chunksize, n_gates])
# import LP model in solver
s = CyClpSimplex(model)
# disable logging in multiprocessing anyway
s.logLevel = 0
i = 0
for raynum in range(ray, ray + chunksize):
# set new B_vector values for actual ray
s.setRowLowerArray(np.squeeze(np.asarray(B_vectors[raynum])))
# set new weights (objectives) for actual ray
s.setObjectiveArray(np.squeeze(np.asarray(weights[raynum])))
# solve with dual method, it is faster
s.dual()
# extract primal solution
soln[i, :] = s.primalVariableSolution["x"][n_gates : 2 * n_gates]
i = i + 1
return soln
def LP_solver_cylp_mp(A_Matrix, B_vectors, weights, really_verbose=False, proc=1):
"""
Solve the Linear Programming problem given in Giangrande et al, 2012 using
the CyLP module using multiple processes.
Parameters
----------
A_Matrix : matrix
Row augmented A matrix, see :py:func:`construct_A_matrix`
B_vectors : matrix
Matrix containing B vectors, see :py:func:`construct_B_vectors`
weights : array
Weights.
really_verbose : bool, optional
True to print CLP messaging. False to suppress.
proc : int, optional
Number of worker processes.
Returns
-------
soln : array
Solution to LP problem.
See Also
--------
LP_solver_cvxopt : Solve LP problem using the CVXOPT module.
LP_solver_pyglpk : Solve LP problem using the PyGLPK module.
LP_solver_cylp : Solve LP problem using the CyLP module using single
process.
"""
import multiprocessing as mp
from cylp.py.modeling.CyLPModel import CyLPArray, CyLPModel
n_gates = weights.shape[1] // 2
n_rays = B_vectors.shape[0]
soln = np.zeros([n_rays, n_gates])
# Create CyLPModel and initialize it
model = CyLPModel()
G = np.matrix(A_Matrix)
h = CyLPArray(np.empty(B_vectors.shape[1]))
x = model.addVariable("x", G.shape[1])
model.addConstraint(G * x >= h)
c = CyLPArray(np.empty(weights.shape[1]))
model.objective = c * x
chunksize = int(n_rays / proc)
# check if equal sized chunks can be distributed to worker processes
if n_rays % chunksize != 0:
print(
"Problem of %d rays cannot be split to %d worker processes!\n\r"
"Fallback to 1 process!" % (n_rays, proc)
)
chunksize = n_rays # fall back to one process
proc = 1
print("Calculating with %d processes, %d rays per chunk" % (proc, chunksize))
def worker(model, B_vectors, weights, ray, chunksize, out_q):
"""
The worker function, invoked in a process.
The results are placed in a dictionary that's pushed to a queue.
"""
outdict = {}
iray = int(ray / chunksize)
outdict[iray] = solve_cylp(model, B_vectors, weights, ray, chunksize)
out_q.put(outdict)
# Queue for LP solutions
out_q = mp.Queue()
procs = []
# fire off worker processes
for raynum in range(0, n_rays, chunksize):
p = mp.Process(
target=worker, args=(model, B_vectors, weights, raynum, chunksize, out_q)
)
procs.append(p)
p.start()
# collecting results
resultdict = {}
for raynum in range(0, n_rays, chunksize):
resultdict.update(out_q.get())
# Wait for all worker processes to finish
for p in procs:
p.join()
# copy results in output array
for raynum in range(0, int(n_rays / chunksize)):
soln[raynum * chunksize : raynum * chunksize + chunksize, :] = resultdict[
raynum
]
# apply smoothing filter to output array
soln = smooth_and_trim_scan(soln, window_len=5, window="sg_smooth")
return soln
def LP_solver_cylp(A_Matrix, B_vectors, weights, really_verbose=False):
"""
Solve the Linear Programming problem given in Giangrande et al, 2012 using
the CyLP module.
Parameters
----------
A_Matrix : matrix
Row augmented A matrix, see :py:func:`construct_A_matrix`
B_vectors : matrix
Matrix containing B vectors, see :py:func:`construct_B_vectors`
weights : array
Weights.
really_verbose : bool, optional
True to print CLP messaging. False to suppress.
Returns
-------
soln : array
Solution to LP problem.
See Also
--------
LP_solver_cvxopt : Solve LP problem using the CVXOPT module.
LP_solver_pyglpk : Solve LP problem using the PyGLPK module.
"""
from cylp.cy.CyClpSimplex import CyClpSimplex
from cylp.py.modeling.CyLPModel import CyLPArray, CyLPModel
n_gates = weights.shape[1] // 2
n_rays = B_vectors.shape[0]
soln = np.zeros([n_rays, n_gates])
# Create CyLPModel and initialize it
model = CyLPModel()
G = np.matrix(A_Matrix)
h = CyLPArray(np.empty(B_vectors.shape[1]))
x = model.addVariable("x", G.shape[1])
model.addConstraint(G * x >= h)
c = CyLPArray(np.squeeze(weights[0]))
model.objective = c * x
# import model in solver
s = CyClpSimplex(model)
# disable logging
if not really_verbose:
s.logLevel = 0
for raynum in range(n_rays):
# set new B_vector values for actual ray
s.setRowLowerArray(np.squeeze(np.asarray(B_vectors[raynum])))
# set new weights (objectives) for actual ray
# solve with dual method, it is faster
s.dual()
# extract primal solution
soln[raynum, :] = s.primalVariableSolution["x"][n_gates : 2 * n_gates]
# apply smoothing filter on a per scan basis
soln = smooth_and_trim_scan(soln, window_len=5, window="sg_smooth")
return soln
[docs]def phase_proc_lp(
radar,
offset,
debug=False,
self_const=60000.0,
low_z=10.0,
high_z=53.0,
min_phidp=0.01,
min_ncp=0.5,
min_rhv=0.8,
fzl=4000.0,
sys_phase=0.0,
overide_sys_phase=False,
nowrap=None,
really_verbose=False,
LP_solver="pyglpk",
refl_field=None,
ncp_field=None,
rhv_field=None,
phidp_field=None,
kdp_field=None,
unf_field=None,
window_len=35,
proc=1,
coef=0.914,
):
"""
Phase process using a LP method [1].
Parameters
----------
radar : Radar
Input radar.
offset : float
Reflectivity offset in dBz.
debug : bool, optional
True to print debugging information.
self_const : float, optional
Self consistency factor.
low_z : float, optional
Low limit for reflectivity. Reflectivity below this value is set to
this limit.
high_z : float, optional
High limit for reflectivity. Reflectivity above this value is set to
this limit.
min_phidp : float, optional
Minimum Phi differential phase.
min_ncp : float, optional
Minimum normal coherent power.
min_rhv : float, optional
Minimum copolar coefficient.
fzl : float, optional
Maximum altitude.
sys_phase : float, optional
System phase in degrees.
overide_sys_phase : bool, optional
True to use `sys_phase` as the system phase. False will calculate a
value automatically.
nowrap : int or None, optional
Gate number to begin phase unwrapping. None will unwrap all phases.
really_verbose : bool, optional
True to print LPX messaging. False to suppress.
LP_solver : 'pyglpk' or 'cvxopt', 'cylp', or 'cylp_mp', optional
Module to use to solve LP problem. Default is 'pyglpk'.
refl_field, ncp_field, rhv_field, phidp_field, kdp_field : str, optional
Name of field in radar which contains the horizonal reflectivity,
normal coherent power, copolar coefficient, differential phase shift,
and differential phase. A value of None for any of these parameters
will use the default field name as defined in the Py-ART configuration
file.
unf_field : str, optional
Name of field which will be added to the radar object which will
contain the unfolded differential phase. Metadata for this field
will be taken from the phidp_field. A value of None will use
the default field name as defined in the Py-ART configuration file.
window_len : int, optional
Length of Sobel window applied to PhiDP field when prior to
calculating KDP.
proc : int, optional
Number of worker processes, only used when `LP_solver` is 'cylp_mp'.
coef : float, optional
Exponent linking Z to KDP in self consistency. kdp=(10**(0.1z))*coef
Returns
-------
reproc_phase : dict
Field dictionary containing processed differential phase shifts.
sob_kdp : dict
Field dictionary containing recalculated differential phases.
References
----------
[1] Giangrande, S.E., R. McGraw, and L. Lei. An Application of
Linear Programming to Polarimetric Radar Differential Phase Processing.
J. Atmos. and Oceanic Tech, 2013, 30, 1716.
"""
# parse the field parameters
if refl_field is None:
refl_field = get_field_name("reflectivity")
if ncp_field is None:
ncp_field = get_field_name("normalized_coherent_power")
if rhv_field is None:
rhv_field = get_field_name("cross_correlation_ratio")
if phidp_field is None:
phidp_field = get_field_name("differential_phase")
if kdp_field is None:
kdp_field = get_field_name("specific_differential_phase")
if unf_field is None:
unf_field = get_field_name("unfolded_differential_phase")
# prepare reflectivity field
refl = copy.deepcopy(radar.fields[refl_field]["data"]) + offset
is_low_z = (refl) < low_z
is_high_z = (refl) > high_z
refl[np.where(is_high_z)] = high_z
refl[np.where(is_low_z)] = low_z
z_mod = refl
# unfold Phi_DP
if debug:
print("Unfolding")
my_unf = get_phidp_unf(
radar,
ncp_lev=min_ncp,
rhohv_lev=min_rhv,
debug=debug,
ncpts=2,
doc=None,
sys_phase=sys_phase,
nowrap=nowrap,
overide_sys_phase=overide_sys_phase,
refl_field=refl_field,
ncp_field=ncp_field,
rhv_field=rhv_field,
phidp_field=phidp_field,
)
my_new_ph = copy.deepcopy(radar.fields[phidp_field])
my_unf[:, -1] = my_unf[:, -2]
my_new_ph["data"] = my_unf
radar.fields.update({unf_field: my_new_ph})
phidp_mod = copy.deepcopy(radar.fields[unf_field]["data"])
phidp_neg = phidp_mod < min_phidp
phidp_mod[np.where(phidp_neg)] = min_phidp
# process
proc_ph = copy.deepcopy(radar.fields[phidp_field])
proc_ph["data"] = phidp_mod
St_Gorlv_differential_5pts = [-0.2, -0.1, 0, 0.1, 0.2]
for sweep in range(len(radar.sweep_start_ray_index["data"])):
if debug:
print("Doing ", sweep)
end_gate, start_ray, end_ray = det_process_range(radar, sweep, fzl, doc=15)
start_gate = 0
A_Matrix = construct_A_matrix(
len(radar.range["data"][start_gate:end_gate]), St_Gorlv_differential_5pts
)
B_vectors = construct_B_vectors(
phidp_mod[start_ray:end_ray, start_gate:end_gate],
z_mod[start_ray:end_ray, start_gate:end_gate],
St_Gorlv_differential_5pts,
dweight=self_const,
coef=coef,
)
weights = np.ones(phidp_mod[start_ray:end_ray, start_gate:end_gate].shape)
nw = np.bmat([weights, np.zeros(weights.shape)])
if LP_solver == "pyglpk":
mysoln = LP_solver_pyglpk(
A_Matrix, B_vectors, nw, really_verbose=really_verbose
)
elif LP_solver == "cvxopt":
mysoln = LP_solver_cvxopt(A_Matrix, B_vectors, nw)
elif LP_solver == "cylp":
mysoln = LP_solver_cylp(
A_Matrix, B_vectors, nw, really_verbose=really_verbose
)
elif LP_solver == "cylp_mp":
mysoln = LP_solver_cylp_mp(
A_Matrix, B_vectors, nw, really_verbose=really_verbose, proc=proc
)
else:
raise ValueError("unknown LP_solver:" + LP_solver)
proc_ph["data"][start_ray:end_ray, start_gate:end_gate] = mysoln
last_gates = proc_ph["data"][start_ray:end_ray, -16]
proc_ph["data"][start_ray:end_ray, -16:] = np.meshgrid(np.ones([16]), last_gates)[1]
proc_ph["valid_min"] = 0.0 # XXX is this correct?
proc_ph["valid_max"] = 400.0 # XXX is this correct?
# prepare output
sobel = 2.0 * np.arange(window_len) / (window_len - 1.0) - 1.0
sobel = sobel / (abs(sobel).sum())
sobel = sobel[::-1]
gate_spacing = (radar.range["data"][1] - radar.range["data"][0]) / 1000.0
kdp = scipy.ndimage.convolve1d(proc_ph["data"], sobel, axis=1) / (
(window_len / 3.0) * 2.0 * gate_spacing
)
# copy the KDP metadata from existing field or create anew
if kdp_field in radar.fields:
sob_kdp = copy.deepcopy(radar.fields[kdp_field])
else:
sob_kdp = get_metadata(kdp_field)
sob_kdp["data"] = kdp
sob_kdp["_FillValue"] = get_fillvalue()
return proc_ph, sob_kdp
[docs]def phase_proc_lp_gf(
radar,
gatefilter=None,
debug=False,
self_const=60000.0,
low_z=10.0,
high_z=53.0,
min_phidp=0.01,
fzl=4000.0,
system_phase=None,
nowrap=None,
really_verbose=False,
LP_solver="pyglpk",
refl_field=None,
phidp_field=None,
kdp_field=None,
unf_field=None,
window_len=35,
proc=1,
coef=0.914,
ncpts=None,
first_gate_sysp=None,
offset=0.0,
doc=0,
):
"""
Phase process using a LP method [1] using Py-ART's Gatefilter.
Parameters
----------
radar : Radar
Input radar.
gatefilter : Gatefilter, optional
Py-ART gatefilter object indicating where processing should be
carried out
debug : bool, optional
True to print debugging information.
self_const : float, optional
Self consistency factor.
low_z : float, optional
Low limit for reflectivity. Reflectivity below this value is set to
this limit.
high_z : float, optional
High limit for reflectivity. Reflectivity above this value is set to
this limit.
fzl : float, optional
Maximum altitude.
system_phase : float, optional
System phase in degrees.
nowrap : int or None, optional
Gate number to begin phase unwrapping. None will unwrap all phases.
really_verbose : bool, optional
True to print LPX messaging. False to suppress.
LP_solver : 'pyglpk' or 'cvxopt', 'cylp', or 'cylp_mp', optional
Module to use to solve LP problem. Default is 'pyglpk'.
refl_field, ncp_field, rhv_field, phidp_field, kdp_field : str, optional
Name of field in radar which contains the horizonal reflectivity,
normal coherent power, copolar coefficient, differential phase shift,
and differential phase. A value of None for any of these parameters
will use the default field name as defined in the Py-ART configuration
file.
unf_field : str, optional
Name of field which will be added to the radar object which will
contain the unfolded differential phase. Metadata for this field
will be taken from the phidp_field. A value of None will use
the default field name as defined in the Py-ART configuration file.
window_len : int, optional
Length of Sobel window applied to PhiDP field when prior to
calculating KDP.
proc : int, optional
Number of worker processes, only used when `LP_solver` is 'cylp_mp'.
coef : float, optional
Exponent linking Z to KDP in self consistency. kdp=(10**(0.1z))*coef
ncpts : int, optional
Minimum number of points in a ray. Regions within a ray smaller than
this or beginning before this gate number are excluded from unfolding.
offset : float, optional
Reflectivity offset to add in dBz.
doc : int, optional
Number of gates to "doc" off the end of a ray.
Returns
-------
reproc_phase : dict
Field dictionary containing processed differential phase shifts.
sob_kdp : dict
Field dictionary containing recalculated differential phases.
References
----------
[1] Giangrande, S.E., R. McGraw, and L. Lei. An Application of
Linear Programming to Polarimetric Radar Differential Phase Processing.
J. Atmos. and Oceanic Tech, 2013, 30, 1716.
"""
# parse the field parameters
if refl_field is None:
refl_field = get_field_name("reflectivity")
if phidp_field is None:
phidp_field = get_field_name("differential_phase")
if kdp_field is None:
kdp_field = get_field_name("specific_differential_phase")
if unf_field is None:
unf_field = get_field_name("unfolded_differential_phase")
# if there is no gatefilter included, include all
if gatefilter is None:
gatefilter = GateFilter(radar)
gatefilter.include_all()
# prepare reflectivity field
refl = copy.deepcopy(radar.fields[refl_field]["data"]) + offset
is_low_z = (refl) < low_z
is_high_z = (refl) > high_z
refl[np.where(is_high_z)] = high_z
refl[np.where(is_low_z)] = low_z
z_mod = refl
# unfold Phi_DP
if debug:
print("Unfolding")
my_unf = get_phidp_unf_gf(
radar,
gatefilter,
debug=debug,
ncpts=ncpts,
sys_phase=system_phase,
nowrap=nowrap,
phidp_field=phidp_field,
first_gate_sysp=first_gate_sysp,
)
my_new_ph = copy.deepcopy(radar.fields[phidp_field])
my_unf[:, -1] = my_unf[:, -2]
my_new_ph["data"] = my_unf
radar.fields.update({unf_field: my_new_ph})
phidp_mod = copy.deepcopy(radar.fields[unf_field]["data"])
phidp_neg = phidp_mod < min_phidp
phidp_mod[np.where(phidp_neg)] = min_phidp
# process
proc_ph = copy.deepcopy(radar.fields[phidp_field])
proc_ph["data"] = phidp_mod
St_Gorlv_differential_5pts = [-0.2, -0.1, 0, 0.1, 0.2]
for sweep in range(len(radar.sweep_start_ray_index["data"])):
if debug:
print("Doing ", sweep)
end_gate, start_ray, end_ray = det_process_range(radar, sweep, fzl, doc=doc)
start_gate = 0
if len(radar.range["data"][start_gate:end_gate]) <= 5:
mysoln = np.zeros_like(
proc_ph["data"][start_ray:end_ray, start_gate:end_gate]
)
proc_ph["data"][start_ray:end_ray, start_gate:end_gate] = mysoln
continue
A_Matrix = construct_A_matrix(
len(radar.range["data"][start_gate:end_gate]), St_Gorlv_differential_5pts
)
B_vectors = construct_B_vectors(
phidp_mod[start_ray:end_ray, start_gate:end_gate],
z_mod[start_ray:end_ray, start_gate:end_gate],
St_Gorlv_differential_5pts,
dweight=self_const,
coef=coef,
)
weights = np.ones(phidp_mod[start_ray:end_ray, start_gate:end_gate].shape)
nw = np.bmat([weights, np.zeros(weights.shape)])
if LP_solver == "pyglpk":
mysoln = LP_solver_pyglpk(
A_Matrix, B_vectors, nw, really_verbose=really_verbose
)
elif LP_solver == "cvxopt":
mysoln = LP_solver_cvxopt(A_Matrix, B_vectors, nw)
elif LP_solver == "cylp":
mysoln = LP_solver_cylp(
A_Matrix, B_vectors, nw, really_verbose=really_verbose
)
elif LP_solver == "cylp_mp":
mysoln = LP_solver_cylp_mp(
A_Matrix, B_vectors, nw, really_verbose=really_verbose, proc=proc
)
else:
raise ValueError("unknown LP_solver:" + LP_solver)
proc_ph["data"][start_ray:end_ray, start_gate:end_gate] = mysoln
last_gates = proc_ph["data"][start_ray:end_ray, -16]
proc_ph["data"][start_ray:end_ray, -16:] = np.meshgrid(np.ones([16]), last_gates)[1]
proc_ph["valid_min"] = 0.0 # XXX is this correct?
proc_ph["valid_max"] = 400.0 # XXX is this correct?
# prepare output
sobel = 2.0 * np.arange(window_len) / (window_len - 1.0) - 1.0
sobel = sobel / (abs(sobel).sum())
sobel = sobel[::-1]
gate_spacing = (radar.range["data"][1] - radar.range["data"][0]) / 1000.0
kdp = scipy.ndimage.convolve1d(proc_ph["data"], sobel, axis=1) / (
(window_len / 3.0) * 2.0 * gate_spacing
)
# copy the KDP metadata from existing field or create anew
if kdp_field in radar.fields:
sob_kdp = copy.deepcopy(radar.fields[kdp_field])
else:
sob_kdp = get_metadata(kdp_field)
sob_kdp["data"] = kdp
sob_kdp["_FillValue"] = get_fillvalue()
return proc_ph, sob_kdp
def get_phidp_unf_gf(
radar,
gatefilter,
debug=False,
ncpts=2,
sys_phase=None,
nowrap=None,
phidp_field=None,
first_gate_sysp=None,
):
"""
Get Unfolded Phi differential phase in areas not gatefiltered.
Parameters
----------
radar : Radar
The input radar.
gatefilter : GateFilter
Only apply on areas incuded in the gatefilter
debug : bool, optioanl
True to print debugging information, False to supress printing.
ncpts : int, optional
Minimum number of points in a ray. Regions within a ray smaller than
this or beginning before this gate number are excluded from
calculations.
sys_phase : float, optional
System phase overide.
nowrap : int or None, optional
Gate number where unwrapping should begin. `None` will unwrap all
gates.
refl_field ncp_field, rhv_field, phidp_field : str, optional
Field names within the radar object which represent the horizontal
reflectivity, normal coherent power, the copolar coefficient, and the
differential phase shift. A value of None for any of these parameters
will use the default field name as defined in the Py-ART
configuration file.
Returns
-------
cordata : array
Unwrapped phi differential phase.
"""
# parse the field parameters
if phidp_field is None:
phidp_field = get_field_name("differential_phase")
t = time()
my_phidp = radar.fields[phidp_field]["data"]
if sys_phase is not None:
system_zero = sys_phase
else:
system_zero = det_sys_phase_gf(
radar, gatefilter, phidp_field=phidp_field, first_gate=first_gate_sysp
)
if system_zero is None:
system_zero = -135
cordata = np.zeros(my_phidp.shape, dtype=float)
all_non_meteo = gatefilter.gate_excluded
for radial in range(my_phidp.shape[0]):
# deterimine non meteo from gatefilter CHANGED
notmeteo = all_non_meteo[radial, :]
x_ma = ma.masked_where(notmeteo, my_phidp[radial, :])
try:
ma.notmasked_contiguous(x_ma)
for slc in ma.notmasked_contiguous(x_ma):
# so trying to get rid of clutter and small things that
# should not add to phidp anyway
if slc.stop - slc.start < ncpts or slc.start < ncpts:
x_ma.mask[slc.start - 1 : slc.stop + 1] = True
c = 0
except TypeError: # non sequence, no valid regions
c = 1 # ie do nothing
x_ma.mask = True
except AttributeError:
c = 1 # also do nothing
x_ma.mask = True
if nowrap is not None:
# Start the unfolding a bit later in order to avoid false
# jumps based on clutter
unwrapped = copy.deepcopy(x_ma)
end_unwrap = unwrap_masked(x_ma[nowrap::], centered=False)
unwrapped[nowrap::] = end_unwrap
else:
unwrapped = unwrap_masked(x_ma, centered=False)
# end so no clutter expected
system_max = (
unwrapped[np.where(np.logical_not(notmeteo))][-10:-1].mean() - system_zero
)
unwrapped_fixed = np.zeros(len(x_ma), dtype=float)
based = unwrapped - system_zero
based[0] = 0.0
notmeteo[0] = False
based[-1] = system_max
notmeteo[-1] = False
unwrapped_fixed[np.where(np.logical_not(based.mask))[0]] = based[
np.where(np.logical_not(based.mask))[0]
]
if len(based[np.where(np.logical_not(based.mask))[0]]) > 11:
unwrapped_fixed[np.where(based.mask)[0]] = np.interp(
np.where(based.mask)[0],
np.where(np.logical_not(based.mask))[0],
smooth_and_trim(based[np.where(np.logical_not(based.mask))[0]]),
)
else:
unwrapped_fixed[np.where(based.mask)[0]] = np.interp(
np.where(based.mask)[0],
np.where(np.logical_not(based.mask))[0],
based[np.where(np.logical_not(based.mask))[0]],
)
if c != 1:
cordata[radial, :] = unwrapped_fixed
else:
cordata[radial, :] = np.zeros(my_phidp.shape[1])
if debug:
print("Exec time: ", time() - t)
return cordata
def det_sys_phase_gf(radar, gatefilter, phidp_field=None, first_gate=30.0):
"""
Determine the system phase.
Parameters
----------
radar : Radar
Radar object for which to determine the system phase.
gatefilter : Gatefilter
Gatefilter object highlighting valid gates.
phidp_field : str, optional
Field name within the radar object which represents
differential phase shift. A value of None will use the default
field name as defined in the Py-ART configuration file.
first_gate : int, optional
Gate index for where to being applying the gatefilter.
Returns
-------
sys_phase : float or None
Estimate of the system phase. None is not estimate can be made.
"""
# parse the field parameters
if phidp_field is None:
phidp_field = get_field_name("differential_phase")
phidp = radar.fields[phidp_field]["data"][:, first_gate:]
last_ray_idx = radar.sweep_end_ray_index["data"][0]
is_meteo = gatefilter.gate_included[:, first_gate:]
return _det_sys_phase_gf(phidp, last_ray_idx, is_meteo)
def _det_sys_phase_gf(phidp, last_ray_idx, radar_meteo):
"""Determine the system phase, see :py:func:`det_sys_phase`."""
good = False
phases = []
for radial in range(last_ray_idx + 1):
meteo = radar_meteo[radial, :]
mpts = np.where(meteo)
if len(mpts[0]) > 25:
good = True
msmth_phidp = smooth_and_trim(phidp[radial, mpts[0]], 9)
phases.append(msmth_phidp[0:25].min())
if not good:
return None
return np.median(phases)