""" Module containing calculations for spectra moments. """
import numpy as np
from ..config import get_metadata
from ..util.hildebrand_sekhon import estimate_noise_hs74
[docs]def spectra_moments(radar):
"""Retrieves the radar moments using a spectra radar object.
Parameter
---------
radar : RadarSpectra
Radar spectra object to use for the calculations.
Returns
-------
fields : dict
Field dictionaries containing moment data.
"""
field_list = {}
times = len(radar.time.values)
rng = len(radar.range.values)
ref = np.zeros((times, rng))
vel = np.zeros((times, rng))
spec_width = np.zeros((times, rng))
skew = np.zeros((times, rng))
kurt = np.zeros((times, rng))
spectra = radar.ds.spectra.values
velocity_bins = radar.velocity_bins.values
d_spectra = np.zeros((times, rng, len(velocity_bins) * 3))
print(spectra.shape)
wavelength = radar.ds.attrs["wavelength"]
for i in range(times):
if i % 100 == 0:
print("Dealiasing %d/%d" % (i, times))
spectra_idx = spectra[i, :, :]
# Get the raw spectra, bins, and wavelength
# Subtract the noise floor from the spectra. Also, return the integration limits
noise_floor, left_limit, right_limit, the_spectra = _get_noise_floor_and_limits(
spectra_idx, avg_window=5
)
d_spec, dealiased_bins = dealias_spectra(
the_spectra, velocity_bins, wavelength, left_limit, right_limit
)
noise_floor = np.tile(noise_floor, (d_spec.shape[1], 1)).T
d_spectra[i, :, :] = np.where(d_spec < noise_floor, np.nan, d_spec)
del the_spectra, d_spec
ref_dict = get_metadata("reflectivity")
ref = _get_reflectivity(d_spectra, dealiased_bins, wavelength)
ref_dict["data"] = ref
field_list["reflectivity"] = ref_dict
vel_dict = get_metadata("velocity")
vel = _get_mean_velocity(d_spectra, dealiased_bins, wavelength, ref)
vel_dict["data"] = vel
field_list["velocity"] = vel_dict
spec_dict = get_metadata("spectrum_width")
spec_width = _get_spectral_width(d_spectra, dealiased_bins, wavelength, ref, vel)
spec_dict["data"] = spec_width
field_list["spectrum_width"] = spec_dict
skew_dict = {"long_name": "skewness", "standard_name": "skewness"}
skew = _get_skewness(d_spectra, dealiased_bins, wavelength, ref, vel, spec_width)
skew_dict["data"] = skew
skew_dict["coordinates"] = "elevation azimuth range"
field_list["skewness"] = skew_dict
kurt_dict = {"long_name": "kurtosis", "standard_name": "kurtosis"}
kurt = _get_kurtosis(d_spectra, dealiased_bins, wavelength, ref, vel, spec_width)
kurt_dict["data"] = kurt
kurt_dict["coordinates"] = "elevation azimuth range"
field_list["kurtosis"] = kurt_dict
return field_list
[docs]def dealias_spectra(the_spectra, vel_bins, wavelength, left_limit, right_limit):
"""Dealias a spectra.
Parameters
----------
the_spectra : array
Spectra field data to dealias.
vel_bins : array
Velocity bin data.
wavelength : float
Spectra radar wavelength.
Returns
-------
new_spectra : array
Dealiased spectra array.
new_bins : array
New velocity bins from dealiased spectra.
"""
# Calculate mean of gate before to decide whether to dealias left
# or right side of spectrum - continuity check!
ref = _get_reflectivity(the_spectra, vel_bins, wavelength)
mean_vel = _get_mean_velocity(the_spectra, vel_bins, wavelength, ref)
mean_vel = np.array(mean_vel)
new_bins = np.concatenate(
[vel_bins - 2 * vel_bins[-1], vel_bins, vel_bins + 2 * vel_bins[-1]]
)
# Expand interval to -3Vn, 3Vn
n_pts = len(vel_bins)
new_spectra = np.nan * np.ones((the_spectra.shape[0], n_pts * 3))
dealiased_already = np.zeros(the_spectra.shape[0])
# Test for aliasing - look for tail on both right interval and left interval
for i in range(the_spectra.shape[0]):
# First test to dealias: peaks on both sides
# Second test: Discontinuity in velocities
if i > 1:
second_vel = mean_vel[i - 1]
else:
second_vel = mean_vel[i]
if np.isfinite(the_spectra[i, 0]) and np.isfinite(the_spectra[i, -1]):
noise_region = np.where(np.isnan(the_spectra[i]))[0]
if second_vel < 0:
right_tail_len = int(the_spectra.shape[1] - noise_region[-1])
new_spectra[i, n_pts - right_tail_len : n_pts] = the_spectra[
i, n_pts - right_tail_len : n_pts
]
new_spectra[i, n_pts : 2 * n_pts - right_tail_len] = the_spectra[
i, 0 : n_pts - right_tail_len
]
else:
left_tail_len = int(noise_region[0])
new_spectra[i, 2 * n_pts : 2 * n_pts + left_tail_len] = the_spectra[
i, 0:left_tail_len
]
new_spectra[i, n_pts + left_tail_len : 2 * n_pts] = the_spectra[
i, left_tail_len:
]
dealiased_already[i] = 1
# Do a second check for continuity
mean_vel = _get_mean_velocity(new_spectra, new_bins, wavelength, ref)
for i in range(new_spectra.shape[0]):
# First test to dealias: peaks on both sides
# Second test: Discontinuity in velocities
if i > 1:
second_vel = mean_vel[i - 1]
else:
second_vel = mean_vel[i]
# Discontinuity = more than 1 Nyquist switch in mean velocity
if abs(mean_vel[i] - second_vel) > vel_bins[-1] and dealiased_already[i] == 0:
noise_region = np.where(np.isnan(the_spectra[i]))[0]
if mean_vel[i] > 0:
right_tail_len = int(the_spectra.shape[1] - noise_region[-1])
new_spectra[i, n_pts - right_tail_len : n_pts] = the_spectra[
i, n_pts - right_tail_len : n_pts
]
new_spectra[i, n_pts : 2 * n_pts - right_tail_len] = the_spectra[
i, 0 : n_pts - right_tail_len
]
else:
left_tail_len = int(noise_region[0])
new_spectra[i, 2 * n_pts : 2 * n_pts + left_tail_len] = the_spectra[
i, 0:left_tail_len
]
new_spectra[i, n_pts + left_tail_len : 2 * n_pts] = the_spectra[
i, left_tail_len:
]
mean_vel[i] = _get_mean_velocity(
new_spectra[i], new_bins, wavelength, ref[i]
)
dealiased_already[i] = 1
elif dealiased_already[i] == 0:
new_spectra[i, n_pts : 2 * n_pts] = the_spectra[i]
return new_spectra, new_bins
def _get_limits_dealiased_spectra(the_spectra):
"""Calculates limits for a dealiased spectra."""
left = np.zeros(the_spectra.shape[0])
right = np.zeros(the_spectra.shape[0])
new_spec = the_spectra.copy()
for i in range(the_spectra.shape[0]):
try:
peak = np.nanargmax(the_spectra[i])
j = peak
while j > 0 and np.isfinite(the_spectra[i, j]):
j = j - 1
left[i] = j
j = peak
while j < the_spectra.shape[1] - 1 and np.isfinite(the_spectra[i, j]):
j = j + 1
right[i] = j
new_spec[i, 0 : int(left[i]) - 1] = np.nan
new_spec[i, int(right[i]) + 1 : -1] = np.nan
except ValueError:
left[i] = np.nan
right[i] = np.nan
return left, right, new_spec
def _get_reflectivity(spectra, bins, wavelength):
"""Calculates reflectivity from a RadarSpectra object."""
spectra_linear = 10 ** (spectra / 10)
radar_constant = 1e18 * wavelength**4 / (0.93 * np.pi**5)
if len(spectra_linear.shape) == 2:
spec_med = radar_constant * (spectra_linear[:, :-1] + spectra_linear[:, 1:]) / 2
diffs = np.tile(np.diff(bins), (spectra.shape[0], 1))
elif len(spectra_linear.shape) == 3:
spec_med = (
radar_constant * (spectra_linear[:, :, :-1] + spectra_linear[:, :, 1:]) / 2
)
diffs = np.tile(np.diff(bins), (spectra.shape[0], spectra.shape[1], 1))
else:
spec_med = radar_constant * (spectra_linear[:-1] + spectra_linear[1:]) / 2
diffs = np.diff(bins)
ref = np.nansum(spec_med * diffs, axis=-1)
return 10 * np.log10(ref)
def _get_mean_velocity(spectra, bins, wavelength, ref):
"""Calculates mean velocity from a RadarSpectra object."""
spectra_linear = 10 ** (spectra / 10)
radar_constant = 1e18 * wavelength**4 / (0.93 * np.pi**5)
ref = 10 ** (ref / 10)
if len(spectra_linear.shape) == 2:
spec_med = radar_constant * (spectra_linear[:, :-1] + spectra_linear[:, 1:]) / 2
diffs = np.tile(np.diff(bins), (spectra.shape[0], 1))
elif len(spectra_linear.shape) == 3:
spec_med = (
radar_constant * (spectra_linear[:, :, :-1] + spectra_linear[:, :, 1:]) / 2
)
diffs = np.tile(np.diff(bins), (spectra.shape[0], spectra.shape[1], 1))
else:
spec_med = radar_constant * (spectra_linear[:-1] + spectra_linear[1:]) / 2
diffs = np.diff(bins)
bins_med = (bins[:-1] + bins[1:]) / 2.0
if len(spectra_linear.shape) == 2:
bins_med = np.tile(bins_med, (spectra.shape[0], 1))
elif len(spectra_linear.shape) == 3:
bins_med = np.tile(bins_med, (spectra.shape[0], spectra.shape[1], 1))
mean_vel = np.nansum(spec_med * bins_med * diffs, axis=-1) / ref
return mean_vel
def _get_spectral_width(spectra, bins, wavelength, ref, mean_vel):
"""Calculates reflectivity from a RadarSpectra object."""
spectra_linear = 10 ** (spectra / 10)
radar_constant = 1e18 * wavelength**4 / (0.93 * np.pi**5)
ref = 10 ** (ref / 10)
if len(spectra_linear.shape) == 2:
spec_med = radar_constant * (spectra_linear[:, :-1] + spectra_linear[:, 1:]) / 2
diffs = np.tile(np.diff(bins), (spectra.shape[0], 1))
mean_vel = np.tile(mean_vel, (len(bins) - 1, 1)).T
elif len(spectra_linear.shape) == 3:
spec_med = (
radar_constant * (spectra_linear[:, :, :-1] + spectra_linear[:, :, 1:]) / 2
)
diffs = np.tile(np.diff(bins), (spectra.shape[0], spectra.shape[1], 1))
mean_vel = np.tile(mean_vel.T, (len(bins) - 1, 1, 1)).T
else:
spec_med = radar_constant * (spectra_linear[:-1] + spectra_linear[1:]) / 2
diffs = np.diff(bins)
diffs = np.tile(np.diff(bins), (spectra.shape[0], spectra.shape[1], 1))
bins_med = (bins[:-1] + bins[1:]) / 2.0
if len(spectra_linear.shape) == 2:
bins_med = np.tile(bins_med, (spectra.shape[0], 1))
elif len(spectra_linear.shape) == 3:
bins_med = np.tile(bins_med, (spectra.shape[0], spectra.shape[1], 1))
spec_wid = np.nansum(spec_med * (bins_med - mean_vel) ** 2 * diffs, axis=-1) / ref
return np.sqrt(spec_wid)
def _get_skewness(spectra, bins, wavelength, ref, mean_vel, spec_width):
"""Calculates skewness from a RadarSpectra object."""
spectra_linear = 10 ** (spectra / 10)
radar_constant = 1e18 * wavelength**4 / (0.93 * np.pi**5)
ref = 10 ** (ref / 10)
if len(spectra_linear.shape) == 2:
spec_med = radar_constant * (spectra_linear[:, :-1] + spectra_linear[:, 1:]) / 2
diffs = np.tile(np.diff(bins), (spectra.shape[0], 1))
mean_vel = np.tile(mean_vel, (len(bins) - 1, 1)).T
elif len(spectra_linear.shape) == 3:
spec_med = (
radar_constant * (spectra_linear[:, :, :-1] + spectra_linear[:, :, 1:]) / 2
)
diffs = np.tile(np.diff(bins), (spectra.shape[0], spectra.shape[1], 1))
mean_vel = np.tile(mean_vel.T, (len(bins) - 1, 1, 1)).T
else:
spec_med = radar_constant * (spectra_linear[:-1] + spectra_linear[1:]) / 2
diffs = np.diff(bins)
diffs = np.tile(np.diff(bins), (spectra.shape[0], spectra.shape[1], 1))
bins_med = (bins[:-1] + bins[1:]) / 2.0
if len(spectra_linear.shape) == 2:
bins_med = np.tile(bins_med, (spectra.shape[0], 1))
elif len(spectra_linear.shape) == 3:
bins_med = np.tile(bins_med, (spectra.shape[0], spectra.shape[1], 1))
skew = np.nansum(spec_med * (bins_med - mean_vel) ** 3 * diffs, axis=-1) / ref
return skew / spec_width**3
def _get_kurtosis(spectra, bins, wavelength, ref, mean_vel, spec_width):
"""Calculates a Kurtosis field from a RadarSpectra object."""
spectra_linear = 10 ** (spectra / 10)
radar_constant = 1e18 * wavelength**4 / (0.93 * np.pi**5)
ref = 10 ** (ref / 10)
if len(spectra_linear.shape) == 2:
spec_med = radar_constant * (spectra_linear[:, :-1] + spectra_linear[:, 1:]) / 2
diffs = np.tile(np.diff(bins), (spectra.shape[0], 1))
mean_vel = np.tile(mean_vel, (len(bins) - 1, 1)).T
elif len(spectra_linear.shape) == 3:
spec_med = (
radar_constant * (spectra_linear[:, :, :-1] + spectra_linear[:, :, 1:]) / 2
)
diffs = np.tile(np.diff(bins), (spectra.shape[0], spectra.shape[1], 1))
mean_vel = np.tile(mean_vel.T, (len(bins) - 1, 1, 1)).T
else:
spec_med = radar_constant * (spectra_linear[:-1] + spectra_linear[1:]) / 2
diffs = np.diff(bins)
bins_med = (bins[:-1] + bins[1:]) / 2.0
if len(spectra_linear.shape) == 2:
bins_med = np.tile(bins_med, (spectra.shape[0], 1))
elif len(spectra_linear.shape) == 3:
bins_med = np.tile(bins_med, (spectra.shape[0], spectra.shape[1], 1))
kurt = np.nansum(spec_med * (bins_med - mean_vel) ** 4 * diffs, axis=-1) / ref
return kurt / spec_width**4
def _get_noise_floor_and_limits(the_spectra, avg_window=1):
"""Calculates noise floor and limits from a RadarSpectra object."""
lin_spectra = 10 ** (the_spectra / 10.0)
if avg_window > 1:
for i in range(lin_spectra.shape[0]):
lin_spectra[i] = np.convolve(
lin_spectra[i], np.ones(avg_window) / avg_window, mode="same"
)
noise_floor_thresh = np.zeros((lin_spectra.shape[0],))
for i in range(lin_spectra.shape[0]):
noise_floor = estimate_noise_hs74(lin_spectra[i], navg=avg_window)
noise_floor_thresh[i] = 10 * np.log10(noise_floor[0])
the_spectra[i] = 10 * np.log10(lin_spectra[i] - noise_floor[0])
the_spectra[lin_spectra < 0] = np.nan
left = np.nan * np.ones(the_spectra.shape[0])
right = np.nan * np.ones(the_spectra.shape[0])
for i in range(the_spectra.shape[0]):
try:
peak_loc = np.nanargmax(the_spectra[i])
j = peak_loc
while np.isfinite(the_spectra[i, j]) and j > 0:
j = j - 1
left[i] = j
j = peak_loc
while np.isfinite(the_spectra[i, j]) and j < the_spectra.shape[1] - 1:
j = j + 1
right[i] = j
except ValueError:
left[i] = np.nan
right[i] = np.nan
spectra = np.ma.masked_where(np.isnan(the_spectra), the_spectra)
return noise_floor_thresh, left, right, spectra