"""
pyart.retrieve.quasi_vertical_profile
=====================================
Retrieval of QVPs from a radar object
.. autosummary::
:toctree: generated/
quasi_vertical_profile
compute_qvp
compute_rqvp
compute_evp
compute_svp
compute_vp
compute_ts_along_coord
project_to_vertical
find_rng_index
get_target_elevations
find_nearest_gate
find_neighbour_gates
_create_qvp_object
_create_along_coord_object
_update_qvp_metadata
_update_along_coord_metadata
""" "" "" "" "" "" "" "" ""
from copy import deepcopy
from warnings import warn
import numpy as np
from netCDF4 import num2date
from scipy.interpolate import interp1d
from ..core.transforms import antenna_to_cartesian
from ..io.common import make_time_unit_str
from ..util.circular_stats import compute_directional_stats
from ..util.datetime_utils import datetime_from_radar
from ..util.radar_utils import ma_broadcast_to
from ..util.xsect import cross_section_ppi, cross_section_rhi
[docs]def quasi_vertical_profile(radar, desired_angle=None, fields=None, gatefilter=None):
"""
Quasi Vertical Profile.
Creates a QVP object containing fields from a radar object that can
be used to plot and produce the quasi vertical profile
Parameters
----------
radar : Radar
Radar object used.
field : string
Radar field to use for QVP calculation.
desired_angle : float
Radar tilt angle to use for indexing radar field data.
None will result in wanted_angle = 20.0
Other Parameters
----------------
gatefilter : GateFilter
A GateFilter indicating radar gates that should be excluded
from the import qvp calculation
Returns
-------
qvp : Dictonary
A quasi vertical profile object containing fields
from a radar object
References
----------
Troemel, S., M. Kumjian, A. Ryzhkov, and C. Simmer, 2013: Backscatter
differential phase - estimation and variability. J Appl. Meteor. Clim.
52, 2529 - 2548.
Troemel, S., A. Ryzhkov, P. Zhang, and C. Simmer, 2014:
Investigations of backscatter differential phase in the melting layer. J.
Appl. Meteorol. Clim. 54, 2344 - 2359.
Ryzhkov, A., P. Zhang, H. Reeves, M. Kumjian, T. Tschallener, S. Tromel,
C. Simmer, 2015: Quasi-vertical profiles - a new way to look at
polarimetric radar data. Submitted to J. Atmos. Oceanic Technol.
"""
# Creating an empty dictonary
qvp = {}
# Setting the desired radar angle and getting index value for desired
# radar angle
if desired_angle is None:
desired_angle = 20.0
index = abs(radar.fixed_angle["data"] - desired_angle).argmin()
radar_slice = radar.get_slice(index)
# Setting field parameters
# If fields is None then all radar fields pulled else defined field is used
if fields is None:
fields = radar.fields
for field in fields:
# Filtering data based on defined gatefilter
# If none is defined goes to else statement
if gatefilter is not None:
get_fields = radar.get_field(index, field)
mask_fields = np.ma.masked_where(
gatefilter.gate_excluded[radar_slice], get_fields
)
radar_fields = np.ma.mean(mask_fields, axis=0)
else:
radar_fields = radar.get_field(index, field).mean(axis=0)
qvp.update({field: radar_fields})
else:
# Filtering data based on defined gatefilter
# If none is defined goes to else statement
if gatefilter is not None:
get_field = radar.get_field(index, fields)
mask_field = np.ma.masked_where(
gatefilter.gate_excluded[radar_slice], get_field
)
radar_field = np.ma.mean(mask_field, axis=0)
else:
radar_field = radar.get_field(index, fields).mean(axis=0)
qvp.update({fields: radar_field})
# Adding range, time, and height fields
qvp.update({"range": radar.range["data"], "time": radar.time})
_, _, z = antenna_to_cartesian(
qvp["range"] / 1000.0, 0.0, radar.fixed_angle["data"][index]
)
qvp.update({"height": z})
return qvp
[docs]def compute_qvp(
radar,
field_names,
ref_time=None,
angle=0.0,
ang_tol=1.0,
hmax=10000.0,
hres=50.0,
avg_type="mean",
nvalid_min=30,
interp_kind="none",
qvp=None,
):
"""
Computes quasi vertical profiles.
Parameters
----------
radar : Radar
Radar object used.
field_names : list of str
list of field names to add to the QVP
ref_time : datetime object
reference time for current radar volume
angle : int or float
If the radar object contains a PPI volume, the sweep number to
use, if it contains an RHI volume the elevation angle.
ang_tol : float
If the radar object contains an RHI volume, the tolerance in the
elevation angle for the conversion into PPI
hmax : float
The maximum height to plot [m].
hres : float
The height resolution [m].
avg_type : str
The type of averaging to perform. Can be either "mean" or "median"
nvalid_min : int
Minimum number of valid points to accept average.
interp_kind : str
type of interpolation when projecting to vertical grid: 'none',
or 'nearest', etc.
'none' will select from all data points within the regular grid
height bin the closest to the center of the bin.
'nearest' will select the closest data point to the center of the
height bin regardless if it is within the height bin or not.
Data points can be masked values
If another type of interpolation is selected masked values will be
eliminated from the data points before the interpolation
qvp : QVP object or None
If it is not None this is the QVP object where to store the data from
the current time step. Otherwise a new QVP object will be created
Returns
-------
qvp : qvp object
The computed QVP object
Reference
---------
Ryzhkov A., Zhang P., Reeves H., Kumjian M., Tschallener T., Trömel S.,
Simmer C. 2016: Quasi-Vertical Profiles: A New Way to Look at Polarimetric
Radar Data. JTECH vol. 33 pp 551-562
"""
if avg_type not in ("mean", "median"):
warn("Unsuported statistics " + avg_type)
return None
radar_aux = deepcopy(radar)
# transform radar into ppi over the required elevation
if radar_aux.scan_type == "rhi":
radar_aux = cross_section_rhi(radar_aux, [angle], el_tol=ang_tol)
elif radar_aux.scan_type == "ppi":
radar_aux = radar_aux.extract_sweeps([int(angle)])
else:
warn("Error: unsupported scan type.")
return None
if qvp is None:
qvp = _create_qvp_object(
radar_aux,
field_names,
qvp_type="qvp",
start_time=ref_time,
hmax=hmax,
hres=hres,
)
# modify metadata
if ref_time is None:
ref_time = datetime_from_radar(radar_aux)
qvp = _update_qvp_metadata(
qvp,
ref_time,
qvp.longitude["data"][0],
qvp.latitude["data"][0],
elev=qvp.fixed_angle["data"][0],
)
for field_name in field_names:
# compute QVP data
if field_name not in radar_aux.fields:
warn("Field " + field_name + " not in radar object")
qvp_data = np.ma.masked_all(qvp.ngates)
else:
values, _ = compute_directional_stats(
radar_aux.fields[field_name]["data"],
avg_type=avg_type,
nvalid_min=nvalid_min,
axis=0,
)
# Project to vertical grid:
qvp_data = project_to_vertical(
values,
radar_aux.gate_altitude["data"][0, :],
qvp.range["data"],
interp_kind=interp_kind,
)
# Put data in radar object
if np.size(qvp.fields[field_name]["data"]) == 0:
qvp.fields[field_name]["data"] = qvp_data.reshape(1, qvp.ngates)
else:
qvp.fields[field_name]["data"] = np.ma.concatenate(
(qvp.fields[field_name]["data"], qvp_data.reshape(1, qvp.ngates))
)
return qvp
[docs]def compute_rqvp(
radar,
field_names,
ref_time=None,
hmax=10000.0,
hres=2.0,
avg_type="mean",
nvalid_min=30,
interp_kind="nearest",
rmax=50000.0,
weight_power=2.0,
qvp=None,
):
"""
Computes range-defined quasi vertical profiles.
Parameters
----------
radar : Radar
Radar object used.
field_names : list of str
list of field names to add to the QVP
ref_time : datetime object
reference time for current radar volume
hmax : float
The maximum height to plot [m].
hres : float
The height resolution [m].
avg_type : str
The type of averaging to perform. Can be either "mean" or "median"
nvalid_min : int
Minimum number of valid points to accept average.
interp_kind : str
type of interpolation when projecting to vertical grid: 'none',
or 'nearest', etc.
'none' will select from all data points within the regular grid
height bin the closest to the center of the bin.
'nearest' will select the closest data point to the center of the
height bin regardless if it is within the height bin or not.
Data points can be masked values
If another type of interpolation is selected masked values will be
eliminated from the data points before the interpolation
rmax : float
ground range up to which the data is intended for use [m].
weight_power : float
Power p of the weighting function 1/abs(grng-(rmax-1))**p given to
the data outside the desired range. -1 will set the weight to 0.
qvp : QVP object or None
If it is not None this is the QVP object where to store the data from
the current time step. Otherwise a new QVP object will be created
Returns
-------
qvp : qvp object
The computed range defined quasi vertical profile
Reference
---------
Tobin D.M., Kumjian M.R. 2017: Polarimetric Radar and Surface-Based
Precipitation-Type Observations of ice Pellet to Freezing Rain
Transitions. Weather and Forecasting vol. 32 pp 2065-2082
"""
if avg_type not in ("mean", "median"):
warn("Unsuported statistics " + avg_type)
return None
radar_aux = deepcopy(radar)
# transform radar into ppi over the required elevation
if radar_aux.scan_type == "rhi":
target_elevations, el_tol = get_target_elevations(radar_aux)
radar_ppi = cross_section_rhi(radar_aux, target_elevations, el_tol=el_tol)
elif radar_aux.scan_type == "ppi":
radar_ppi = radar_aux
else:
warn("Error: unsupported scan type.")
return None
if qvp is None:
qvp = _create_qvp_object(
radar_aux,
field_names,
qvp_type="rqvp",
start_time=ref_time,
hmax=hmax,
hres=hres,
)
# modify metadata
if ref_time is None:
ref_time = datetime_from_radar(radar_ppi)
qvp = _update_qvp_metadata(
qvp, ref_time, qvp.longitude["data"][0], qvp.latitude["data"][0], elev=90.0
)
rmax_km = rmax / 1000.0
grng_interp = np.ma.masked_all((radar_ppi.nsweeps, qvp.ngates))
val_interp = dict()
for field_name in field_names:
val_interp.update(
{field_name: np.ma.masked_all((radar_ppi.nsweeps, qvp.ngates))}
)
for sweep in range(radar_ppi.nsweeps):
radar_aux = deepcopy(radar_ppi)
radar_aux = radar_aux.extract_sweeps([sweep])
height = radar_aux.gate_altitude["data"][0, :]
# compute ground range [Km]
grng = (
np.sqrt(
np.power(radar_aux.gate_x["data"][0, :], 2.0)
+ np.power(radar_aux.gate_y["data"][0, :], 2.0)
)
/ 1000.0
)
# Project ground range to grid
f = interp1d(
height, grng, kind=interp_kind, bounds_error=False, fill_value="extrapolate"
)
grng_interp[sweep, :] = f(qvp.range["data"])
for field_name in field_names:
if field_name not in radar_aux.fields:
warn("Field " + field_name + " not in radar object")
continue
# Compute QVP for this sweep
values, _ = compute_directional_stats(
radar_aux.fields[field_name]["data"],
avg_type=avg_type,
nvalid_min=nvalid_min,
axis=0,
)
# Project to grid
val_interp[field_name][sweep, :] = project_to_vertical(
values, height, qvp.range["data"], interp_kind=interp_kind
)
# Compute weight
weight = np.ma.abs(grng_interp - (rmax_km - 1.0))
weight[grng_interp <= rmax_km - 1.0] = 1.0 / np.power(
weight[grng_interp <= rmax_km - 1.0], 0.0
)
if weight_power == -1:
weight[grng_interp > rmax_km - 1.0] = 0.0
else:
weight[grng_interp > rmax_km - 1.0] = 1.0 / np.power(
weight[grng_interp > rmax_km - 1.0], weight_power
)
for field_name in field_names:
# mask weights where there is no data
mask = np.ma.getmaskarray(val_interp[field_name])
weight_aux = np.ma.masked_where(mask, weight)
# Weighted average
qvp_data = np.ma.sum(val_interp[field_name] * weight_aux, axis=0) / np.ma.sum(
weight_aux, axis=0
)
# Put data in radar object
if np.size(qvp.fields[field_name]["data"]) == 0:
qvp.fields[field_name]["data"] = qvp_data.reshape(1, qvp.ngates)
else:
qvp.fields[field_name]["data"] = np.ma.concatenate(
(qvp.fields[field_name]["data"], qvp_data.reshape(1, qvp.ngates))
)
return qvp
[docs]def compute_evp(
radar,
field_names,
lon,
lat,
ref_time=None,
latlon_tol=0.0005,
delta_rng=15000.0,
delta_azi=10,
hmax=10000.0,
hres=250.0,
avg_type="mean",
nvalid_min=1,
interp_kind="none",
qvp=None,
):
"""
Computes enhanced vertical profiles.
Parameters
----------
radar : Radar
Radar object used.
field_names : list of str
list of field names to add to the QVP
lat, lon : float
latitude and longitude of the point of interest [deg]
ref_time : datetime object
reference time for current radar volume
latlon_tol : float
tolerance in latitude and longitude in deg.
delta_rng, delta_azi : float
maximum range distance [m] and azimuth distance [degree] from the
central point of the evp containing data to average.
hmax : float
The maximum height to plot [m].
hres : float
The height resolution [m].
avg_type : str
The type of averaging to perform. Can be either "mean" or "median"
nvalid_min : int
Minimum number of valid points to accept average.
interp_kind : str
type of interpolation when projecting to vertical grid: 'none',
or 'nearest', etc.
'none' will select from all data points within the regular grid
height bin the closest to the center of the bin.
'nearest' will select the closest data point to the center of the
height bin regardless if it is within the height bin or not.
Data points can be masked values
If another type of interpolation is selected masked values will be
eliminated from the data points before the interpolation
qvp : QVP object or None
If it is not None this is the QVP object where to store the data from
the current time step. Otherwise a new QVP object will be created
Returns
-------
qvp : qvp object
The computed enhanced vertical profile
Reference
---------
Kaltenboeck R., Ryzhkov A. 2016: A freezing rain storm explored with a
C-band polarimetric weather radar using the QVP methodology.
Meteorologische Zeitschrift vol. 26 pp 207-222
"""
if avg_type not in ("mean", "median"):
warn("Unsuported statistics " + avg_type)
return None
radar_aux = deepcopy(radar)
# transform radar into ppi over the required elevation
if radar_aux.scan_type == "rhi":
target_elevations, el_tol = get_target_elevations(radar_aux)
radar_ppi = cross_section_rhi(radar_aux, target_elevations, el_tol=el_tol)
elif radar_aux.scan_type == "ppi":
radar_ppi = radar_aux
else:
warn("Error: unsupported scan type.")
return None
radar_aux = radar_ppi.extract_sweeps([0])
if qvp is None:
qvp = _create_qvp_object(
radar_aux,
field_names,
qvp_type="evp",
start_time=ref_time,
hmax=hmax,
hres=hres,
)
# modify metadata
if ref_time is None:
ref_time = datetime_from_radar(radar_aux)
qvp = _update_qvp_metadata(qvp, ref_time, lon, lat, elev=90.0)
height = dict()
values = dict()
for field_name in field_names:
height.update({field_name: np.array([], dtype=float)})
values.update({field_name: np.ma.array([], dtype=float)})
for sweep in range(radar_ppi.nsweeps):
radar_aux = deepcopy(radar_ppi)
radar_aux = radar_aux.extract_sweeps([sweep])
# find nearest gate to lat lon point
ind_ray, _, azi, rng = find_nearest_gate(
radar_aux, lat, lon, latlon_tol=latlon_tol
)
if ind_ray is None:
continue
# find neighbouring gates to be selected
inds_ray, inds_rng = find_neighbour_gates(
radar_aux, azi, rng, delta_azi=delta_azi, delta_rng=delta_rng
)
for field_name in field_names:
if field_name not in radar_aux.fields:
warn("Field " + field_name + " not in radar object")
continue
height[field_name] = np.append(
height[field_name], radar_aux.gate_altitude["data"][ind_ray, inds_rng]
)
# keep only data we are interested in
field = radar_aux.fields[field_name]["data"][:, inds_rng]
field = field[inds_ray, :]
vals, _ = compute_directional_stats(
field, avg_type=avg_type, nvalid_min=nvalid_min, axis=0
)
values[field_name] = np.ma.append(values[field_name], vals)
for field_name in field_names:
# Project to vertical grid:
qvp_data = project_to_vertical(
values[field_name],
height[field_name],
qvp.range["data"],
interp_kind=interp_kind,
)
# Put data in radar object
if np.size(qvp.fields[field_name]["data"]) == 0:
qvp.fields[field_name]["data"] = qvp_data.reshape(1, qvp.ngates)
else:
qvp.fields[field_name]["data"] = np.ma.concatenate(
(qvp.fields[field_name]["data"], qvp_data.reshape(1, qvp.ngates))
)
return qvp
[docs]def compute_svp(
radar,
field_names,
lon,
lat,
angle,
ref_time=None,
ang_tol=1.0,
latlon_tol=0.0005,
delta_rng=15000.0,
delta_azi=10,
hmax=10000.0,
hres=250.0,
avg_type="mean",
nvalid_min=1,
interp_kind="none",
qvp=None,
):
"""
Computes slanted vertical profiles.
Parameters
----------
radar : Radar
Radar object used.
field_names : list of str
list of field names to add to the QVP
lat, lon : float
latitude and longitude of the point of interest [deg]
angle : int or float
If the radar object contains a PPI volume, the sweep number to
use, if it contains an RHI volume the elevation angle.
ref_time : datetime object
reference time for current radar volume
ang_tol : float
If the radar object contains an RHI volume, the tolerance in the
elevation angle for the conversion into PPI
latlon_tol : float
tolerance in latitude and longitude in deg.
delta_rng, delta_azi : float
maximum range distance [m] and azimuth distance [degree] from the
central point of the evp containing data to average.
hmax : float
The maximum height to plot [m].
hres : float
The height resolution [m].
avg_type : str
The type of averaging to perform. Can be either "mean" or "median"
nvalid_min : int
Minimum number of valid points to accept average.
interp_kind : str
type of interpolation when projecting to vertical grid: 'none',
or 'nearest', etc.
'none' will select from all data points within the regular grid
height bin the closest to the center of the bin.
'nearest' will select the closest data point to the center of the
height bin regardless if it is within the height bin or not.
Data points can be masked values
If another type of interpolation is selected masked values will be
eliminated from the data points before the interpolation
qvp : QVP object or None
If it is not None this is the QVP object where to store the data from
the current time step. Otherwise a new QVP object will be created
Returns
-------
qvp : qvp object
The computed slanted vertical profile
Reference
---------
Bukovcic P., Zrnic D., Zhang G. 2017: Winter Precipitation Liquid-Ice
Phase Transitions Revealed with Polarimetric Radar and 2DVD Observations
in Central Oklahoma. JTECH vol. 56 pp 1345-1363
"""
if avg_type not in ("mean", "median"):
warn("Unsuported statistics " + avg_type)
return None
radar_aux = deepcopy(radar)
# transform radar into ppi over the required elevation
if radar_aux.scan_type == "rhi":
radar_aux = cross_section_rhi(radar_aux, [angle], el_tol=ang_tol)
elif radar_aux.scan_type == "ppi":
radar_aux = radar_aux.extract_sweeps([int(angle)])
else:
warn("Error: unsupported scan type.")
return None
if qvp is None:
qvp = _create_qvp_object(
radar_aux,
field_names,
qvp_type="svp",
start_time=ref_time,
hmax=hmax,
hres=hres,
)
# modify metadata
if ref_time is None:
ref_time = datetime_from_radar(radar_aux)
qvp = _update_qvp_metadata(qvp, ref_time, lon, lat, elev=qvp.fixed_angle["data"][0])
# find nearest gate to lat lon point
ind_ray, _, azi, rng = find_nearest_gate(radar_aux, lat, lon, latlon_tol=latlon_tol)
if ind_ray is None:
warn("No data in selected point")
qvp_data = np.ma.masked_all(qvp.ngates)
for field_name in field_names:
# Put data in radar object
if np.size(qvp.fields[field_name]["data"]) == 0:
qvp.fields[field_name]["data"] = qvp_data.reshape(1, qvp.ngates)
else:
qvp.fields[field_name]["data"] = np.ma.concatenate(
(qvp.fields[field_name]["data"], qvp_data.reshape(1, qvp.ngates))
)
return qvp
# find neighbouring gates to be selected
inds_ray, inds_rng = find_neighbour_gates(
radar_aux, azi, rng, delta_azi=delta_azi, delta_rng=delta_rng
)
height = radar_aux.gate_altitude["data"][ind_ray, inds_rng]
for field_name in field_names:
if field_name not in radar_aux.fields:
warn("Field " + field_name + " not in radar object")
qvp_data = np.ma.masked_all(qvp.ngates)
else:
# keep only data we are interested in
field = radar_aux.fields[field_name]["data"][:, inds_rng]
field = field[inds_ray, :]
# compute values
values, _ = compute_directional_stats(
field, avg_type=avg_type, nvalid_min=nvalid_min, axis=0
)
# Project to vertical grid:
qvp_data = project_to_vertical(
values, height, qvp.range["data"], interp_kind=interp_kind
)
# Put data in radar object
if np.size(qvp.fields[field_name]["data"]) == 0:
qvp.fields[field_name]["data"] = qvp_data.reshape(1, qvp.ngates)
else:
qvp.fields[field_name]["data"] = np.ma.concatenate(
(qvp.fields[field_name]["data"], qvp_data.reshape(1, qvp.ngates))
)
return qvp
[docs]def compute_vp(
radar,
field_names,
lon,
lat,
ref_time=None,
latlon_tol=0.0005,
hmax=10000.0,
hres=50.0,
interp_kind="none",
qvp=None,
):
"""
Computes vertical profiles.
Parameters
----------
radar : Radar
Radar object used.
field_names : list of str
list of field names to add to the QVP
lat, lon : float
latitude and longitude of the point of interest [deg]
ref_time : datetime object
reference time for current radar volume
latlon_tol : float
tolerance in latitude and longitude in deg.
hmax : float
The maximum height to plot [m].
hres : float
The height resolution [m].
interp_kind : str
type of interpolation when projecting to vertical grid: 'none',
or 'nearest', etc.
'none' will select from all data points within the regular grid
height bin the closest to the center of the bin.
'nearest' will select the closest data point to the center of the
height bin regardless if it is within the height bin or not.
Data points can be masked values
If another type of interpolation is selected masked values will be
eliminated from the data points before the interpolation
qvp : QVP object or None
If it is not None this is the QVP object where to store the data from
the current time step. Otherwise a new QVP object will be created
Returns
-------
qvp : qvp object
The computed vertical profile
"""
radar_aux = deepcopy(radar)
# transform radar into ppi over the required elevation
if radar_aux.scan_type == "rhi":
target_elevations, el_tol = get_target_elevations(radar_aux)
radar_ppi = cross_section_rhi(radar_aux, target_elevations, el_tol=el_tol)
elif radar_aux.scan_type == "ppi":
radar_ppi = radar_aux
else:
warn("Error: unsupported scan type.")
return None
if qvp is None:
qvp = _create_qvp_object(
radar_ppi,
field_names,
qvp_type="vp",
start_time=ref_time,
hmax=hmax,
hres=hres,
)
# modify metadata
if ref_time is None:
ref_time = datetime_from_radar(radar_ppi)
qvp = _update_qvp_metadata(qvp, ref_time, lon, lat, elev=90.0)
height = dict()
values = dict()
for field_name in field_names:
height.update({field_name: np.array([], dtype=float)})
values.update({field_name: np.ma.array([], dtype=float)})
for sweep in range(radar_ppi.nsweeps):
radar_aux = deepcopy(radar_ppi.extract_sweeps([sweep]))
# find nearest gate to lat lon point
ind_ray, ind_rng, _, _ = find_nearest_gate(
radar_aux, lat, lon, latlon_tol=latlon_tol
)
if ind_ray is None:
continue
for field_name in field_names:
if field_name not in radar_aux.fields:
warn("Field " + field_name + " not in radar object")
continue
height[field_name] = np.append(
height[field_name], radar_aux.gate_altitude["data"][ind_ray, ind_rng]
)
values[field_name] = np.ma.append(
values[field_name],
radar_aux.fields[field_name]["data"][ind_ray, ind_rng],
)
for field_name in field_names:
# Project to vertical grid:
qvp_data = project_to_vertical(
values[field_name],
height[field_name],
qvp.range["data"],
interp_kind=interp_kind,
)
# Put data in radar object
if np.size(qvp.fields[field_name]["data"]) == 0:
qvp.fields[field_name]["data"] = qvp_data.reshape(1, qvp.ngates)
else:
qvp.fields[field_name]["data"] = np.ma.concatenate(
(qvp.fields[field_name]["data"], qvp_data.reshape(1, qvp.ngates))
)
return qvp
[docs]def compute_ts_along_coord(
radar,
field_names,
mode="ALONG_AZI",
fixed_range=None,
fixed_azimuth=None,
fixed_elevation=None,
ang_tol=1.0,
rng_tol=50.0,
value_start=None,
value_stop=None,
ref_time=None,
acoord=None,
):
"""
Computes time series along a particular antenna coordinate, i.e. along
azimuth, elevation or range
Parameters
----------
radar : Radar
Radar object used.
field_names : str
list of field names
mode : str
coordinate to extract data along. Can be ALONG_AZI, ALONG_ELE or
ALONG_RNG
fixed_range, fixed_azimuth, fixed_elevation : float
The fixed range [m], azimuth [deg] or elevation [deg] to extract.
In each mode two of these parameters have to be defined. If they are
not defined they default to 0.
ang_tol, rng_tol : float
The angle tolerance [deg] and range tolerance [m] around the fixed
range or azimuth/elevation
value_start, value_stop : float
The minimum and maximum value at which the data along a coordinate
start and stop
ref_time : datetime object
reference time for current radar volume
acoord : acoord object or None
If it is not None this is the object where to store the data from the
current time step. Otherwise a new acoord object will be created
Returns
-------
acoord : acoord object
The computed data along a coordinate
"""
if mode == "ALONG_RNG":
if value_start is None:
value_start = 0.0
if value_stop is None:
value_stop = radar.range["data"][-1]
if fixed_azimuth is None:
fixed_azimuth = 0.0
if fixed_elevation is None:
fixed_elevation = 0.0
elif mode == "ALONG_AZI":
if value_start is None:
value_start = np.min(radar.azimuth["data"])
if value_stop is None:
value_stop = np.max(radar.azimuth["data"])
if fixed_range is None:
fixed_range = 0.0
if fixed_elevation is None:
fixed_elevation = 0.0
elif mode == "ALONG_ELE":
if value_start is None:
value_start = np.min(radar.elevation["data"])
if value_stop is None:
value_stop = np.max(radar.elevation["data"])
if fixed_range is None:
fixed_range = 0.0
if fixed_azimuth is None:
fixed_azimuth = 0.0
else:
warn("Unknown time series of coordinate mode " + mode)
return None
if mode == "ALONG_RNG":
# rng_values : range
# fixed_angle: elevation
# elevation: elevation
# azimuth: azimuth
vals_dict = {}
for field_name in field_names:
rng_values, vals, _, _ = get_data_along_rng(
radar,
field_name,
[fixed_elevation],
[fixed_azimuth],
ang_tol=ang_tol,
rmin=value_start,
rmax=value_stop,
)
if vals.size == 0:
warn("No data found")
return None
vals_dict.update({field_name: vals})
fixed_angle = fixed_elevation
elevation = fixed_elevation
azimuth = fixed_azimuth
atol = rng_tol
elif mode == "ALONG_AZI":
# rng_values : azimuth
# fixed_angle : elevation
# elevation : elevation
# azimuth : range
vals_dict = {}
for field_name in field_names:
rng_values, vals, _, _ = get_data_along_azi(
radar,
field_name,
[fixed_range],
[fixed_elevation],
rng_tol=rng_tol,
ang_tol=ang_tol,
azi_start=value_start,
azi_stop=value_stop,
)
if vals.size == 0:
warn("No data found")
return None
vals_dict.update({field_name: vals})
fixed_angle = fixed_elevation
elevation = fixed_elevation
azimuth = fixed_range
atol = ang_tol
else:
# rng_values : elevation
# fixed_angle : azimuth
# elevation : range
# azimuth : azimuth
vals_dict = {}
for field_name in field_names:
rng_values, vals, _, _ = get_data_along_ele(
radar,
field_name,
[fixed_range],
[fixed_azimuth],
rng_tol=rng_tol,
ang_tol=ang_tol,
ele_min=value_start,
ele_max=value_stop,
)
if vals.size == 0:
warn("No data found")
return None
vals_dict.update({field_name: vals})
fixed_angle = fixed_azimuth
elevation = fixed_range
azimuth = fixed_azimuth
atol = ang_tol
if acoord is None:
acoord = _create_along_coord_object(
radar, field_names, rng_values, fixed_angle, mode, start_time=ref_time
)
if not np.allclose(rng_values, acoord.range["data"], rtol=1e5, atol=atol):
warn("Unable to add data. xvalues different from previous ones")
return None
# modify metadata
if ref_time is None:
ref_time = datetime_from_radar(radar)
acoord = _update_along_coord_metadata(acoord, ref_time, elevation, azimuth)
# Put data in radar object
for field_name in field_names:
if np.size(acoord.fields[field_name]["data"]) == 0:
acoord.fields[field_name]["data"] = vals_dict[field_name].reshape(
1, acoord.ngates
)
else:
acoord.fields[field_name]["data"] = np.ma.concatenate(
(
acoord.fields[field_name]["data"],
vals_dict[field_name].reshape(1, acoord.ngates),
)
)
return acoord
def project_to_vertical(
data_in, data_height, grid_height, interp_kind="none", fill_value=-9999.0
):
"""
Projects radar data to a regular vertical grid
Parameters
----------
data_in : ndarray 1D
the radar data to project
data_height : ndarray 1D
the height of each radar point
grid_height : ndarray 1D
the regular vertical grid to project to
interp_kind : str
The type of interpolation to use: 'none' or 'nearest'
fill_value : float
The fill value used for interpolation
Returns
-------
data_out : ndarray 1D
The projected data
"""
if data_in.size == 0:
data_out = np.ma.masked_all(grid_height.size)
return data_out
if interp_kind == "none":
hres = grid_height[1] - grid_height[0]
data_out = np.ma.masked_all(grid_height.size)
for ind_r, h in enumerate(grid_height):
ind_h = find_rng_index(data_height, h, rng_tol=hres / 2.0)
if ind_h is None:
continue
data_out[ind_r] = data_in[ind_h]
elif interp_kind == "nearest":
data_filled = data_in.filled(fill_value=fill_value)
f = interp1d(
data_height,
data_filled,
kind=interp_kind,
bounds_error=False,
fill_value=fill_value,
)
data_out = np.ma.masked_values(f(grid_height), fill_value)
else:
valid = np.logical_not(np.ma.getmaskarray(data_in))
height_valid = data_height[valid]
data_valid = data_in[valid]
f = interp1d(
height_valid,
data_valid,
kind=interp_kind,
bounds_error=False,
fill_value=fill_value,
)
data_out = np.ma.masked_values(f(grid_height), fill_value)
return data_out
def find_rng_index(rng_vec, rng, rng_tol=0.0):
"""
Find the range index corresponding to a particular range
Parameters
----------
rng_vec : float array
The range data array where to look for
rng : float
The range to search
rng_tol : float
Tolerance [m]
Returns
-------
ind_rng : int
The range index
"""
dist = np.abs(rng_vec - rng)
ind_rng = np.argmin(dist)
if dist[ind_rng] > rng_tol:
return None
return ind_rng
def get_target_elevations(radar_in):
"""
Gets RHI target elevations
Parameters
----------
radar_in : Radar object
current radar object
Returns
-------
target_elevations : 1D-array
Azimuth angles
el_tol : float
azimuth tolerance
"""
sweep_start = radar_in.sweep_start_ray_index["data"][0]
sweep_end = radar_in.sweep_end_ray_index["data"][0]
target_elevations = np.sort(radar_in.elevation["data"][sweep_start : sweep_end + 1])
el_tol = np.median(target_elevations[1:] - target_elevations[:-1])
return target_elevations, el_tol
def find_nearest_gate(radar, lat, lon, latlon_tol=0.0005):
"""
Find the radar gate closest to a lat,lon point
Parameters
----------
radar : radar object
the radar object
lat, lon : float
The position of the point
latlon_tol : float
The tolerance around this point
Returns
-------
ind_ray, ind_rng : int
The ray and range index
azi, rng : float
the range and azimuth position of the gate
"""
# find gates close to lat lon point
inds_ray_aux, inds_rng_aux = np.where(
np.logical_and(
np.logical_and(
radar.gate_latitude["data"] < lat + latlon_tol,
radar.gate_latitude["data"] > lat - latlon_tol,
),
np.logical_and(
radar.gate_longitude["data"] < lon + latlon_tol,
radar.gate_longitude["data"] > lon - latlon_tol,
),
)
)
if inds_ray_aux.size == 0:
warn(
"No data found at point lat "
+ str(lat)
+ " +- "
+ str(latlon_tol)
+ " lon "
+ str(lon)
+ " +- "
+ str(latlon_tol)
+ " deg"
)
return None, None, None, None
# find closest latitude
ind_min = np.argmin(
np.abs(radar.gate_latitude["data"][inds_ray_aux, inds_rng_aux] - lat)
)
ind_ray = inds_ray_aux[ind_min]
ind_rng = inds_rng_aux[ind_min]
azi = radar.azimuth["data"][ind_ray]
rng = radar.range["data"][ind_rng]
return ind_ray, ind_rng, azi, rng
def find_neighbour_gates(radar, azi, rng, delta_azi=None, delta_rng=None):
"""
Find the neighbouring gates within +-delta_azi and +-delta_rng
Parameters
----------
radar : radar object
the radar object
azi, rng : float
The azimuth [deg] and range [m] of the central gate
delta_azi, delta_rng : float
The extend where to look for
Returns
-------
inds_ray_aux, ind_rng_aux : int
The indices (ray, rng) of the neighbouring gates
"""
# find gates close to lat lon point
if delta_azi is None:
inds_ray = np.ma.arange(radar.azimuth["data"].size)
else:
azi_max = azi + delta_azi
azi_min = azi - delta_azi
if azi_max > 360.0:
azi_max -= 360.0
if azi_min < 0.0:
azi_min += 360.0
if azi_max > azi_min:
inds_ray = np.where(
np.logical_and(
radar.azimuth["data"] < azi_max, radar.azimuth["data"] > azi_min
)
)[0]
else:
inds_ray = np.where(
np.logical_or(
radar.azimuth["data"] > azi_min, radar.azimuth["data"] < azi_max
)
)[0]
if delta_rng is None:
inds_rng = np.ma.arange(radar.range["data"].size)
else:
inds_rng = np.where(
np.logical_and(
radar.range["data"] < rng + delta_rng,
radar.range["data"] > rng - delta_rng,
)
)[0]
return inds_ray, inds_rng
def get_data_along_rng(
radar, field_name, fix_elevations, fix_azimuths, ang_tol=1.0, rmin=None, rmax=None
):
"""
Get data at particular (azimuths, elevations)
Parameters
----------
radar : radar object
the radar object where the data is
field_name : str
name of the field to filter
fix_elevations, fix_azimuths: list of floats
List of elevations, azimuths couples [deg]
ang_tol : float
Tolerance between the nominal angle and the radar angle [deg]
rmin, rmax: float
Min and Max range of the obtained data [m]
Returns
-------
xvals : 1D float array
The ranges of each azi, ele pair
yvals : 1D float array
The values
valid_azi, valid_ele : float arrays
The azi, ele pairs
"""
if rmin is None:
rmin = 0.0
if rmax is None:
rmax = np.max(radar.range["data"])
rng_mask = np.logical_and(radar.range["data"] >= rmin, radar.range["data"] <= rmax)
x = radar.range["data"][rng_mask]
xvals = []
yvals = []
valid_azi = []
valid_ele = []
if radar.scan_type in ("ppi", "vertical_pointing"):
for ele, azi in zip(fix_elevations, fix_azimuths):
if radar.scan_type == "vertical_pointing":
dataset_line = deepcopy(radar)
xvals.extend(x)
else:
ind_sweep = find_ang_index(
radar.fixed_angle["data"], ele, ang_tol=ang_tol
)
if ind_sweep is None:
warn("No elevation angle found for fix_elevation " + str(ele))
continue
new_dataset = radar.extract_sweeps([ind_sweep])
try:
dataset_line = cross_section_ppi(new_dataset, [azi], az_tol=ang_tol)
except OSError:
warn(
" No data found at azimuth "
+ str(azi)
+ " and elevation "
+ str(ele)
)
continue
xvals.append(x)
yvals.append(dataset_line.fields[field_name]["data"][0, rng_mask])
valid_azi.append(dataset_line.azimuth["data"][0])
valid_ele.append(dataset_line.elevation["data"][0])
else:
for ele, azi in zip(fix_elevations, fix_azimuths):
ind_sweep = find_ang_index(radar.fixed_angle["data"], azi, ang_tol=ang_tol)
if ind_sweep is None:
warn("No azimuth angle found for fix_azimuth " + str(azi))
continue
new_dataset = radar.extract_sweeps([ind_sweep])
try:
dataset_line = cross_section_rhi(new_dataset, [ele], el_tol=ang_tol)
except OSError:
warn(
" No data found at azimuth "
+ str(azi)
+ " and elevation "
+ str(ele)
)
continue
yvals.extend(dataset_line.fields[field_name]["data"][0, rng_mask])
xvals.extend(x)
valid_azi.append(dataset_line.azimuth["data"][0])
valid_ele.append(dataset_line.elevation["data"][0])
return np.array(xvals), np.ma.array(yvals), valid_azi, valid_ele
def get_data_along_azi(
radar,
field_name,
fix_ranges,
fix_elevations,
rng_tol=50.0,
ang_tol=1.0,
azi_start=None,
azi_stop=None,
):
"""
Get data at particular (ranges, elevations)
Parameters
----------
radar : radar object
the radar object where the data is
field_name : str
name of the field to filter
fix_ranges, fix_elevations: list of floats
List of ranges [m], elevations [deg] couples
rng_tol : float
Tolerance between the nominal range and the radar range [m]
ang_tol : float
Tolerance between the nominal angle and the radar angle [deg]
azi_start, azi_stop: float
Start and stop azimuth angle of the data [deg]
Returns
-------
xvals : 1D float array
The ranges of each rng, ele pair
yvals : 1D float array
The values
valid_rng, valid_ele : float arrays
The rng, ele pairs
"""
if azi_start is None:
azi_start = np.min(radar.azimuth["data"])
if azi_stop is None:
azi_stop = np.max(radar.azimuth["data"])
yvals = []
xvals = []
valid_rng = []
valid_ele = []
for rng, ele in zip(fix_ranges, fix_elevations):
ind_rng = find_rng_index(radar.range["data"], rng, rng_tol=rng_tol)
if ind_rng is None:
warn("No range gate found for fix_range " + str(rng))
continue
if radar.scan_type == "ppi":
ind_sweep = find_ang_index(radar.fixed_angle["data"], ele, ang_tol=ang_tol)
if ind_sweep is None:
warn("No elevation angle found for fix_elevation " + str(ele))
continue
new_dataset = radar.extract_sweeps([ind_sweep])
else:
try:
new_dataset = cross_section_rhi(radar, [ele], el_tol=ang_tol)
except OSError:
warn(
" No data found at range " + str(rng) + " and elevation " + str(ele)
)
continue
if azi_start < azi_stop:
azi_mask = np.logical_and(
new_dataset.azimuth["data"] >= azi_start,
new_dataset.azimuth["data"] <= azi_stop,
)
else:
azi_mask = np.logical_or(
new_dataset.azimuth["data"] >= azi_start,
new_dataset.azimuth["data"] <= azi_stop,
)
yvals.extend(new_dataset.fields[field_name]["data"][azi_mask, ind_rng])
xvals.extend(new_dataset.azimuth["data"][azi_mask])
valid_rng.append(new_dataset.range["data"][ind_rng])
valid_ele.append(new_dataset.elevation["data"][0])
return np.array(xvals), np.ma.array(yvals), valid_rng, valid_ele
def get_data_along_ele(
radar,
field_name,
fix_ranges,
fix_azimuths,
rng_tol=50.0,
ang_tol=1.0,
ele_min=None,
ele_max=None,
):
"""
Get data at particular (ranges, azimuths)
Parameters
----------
radar : radar object
the radar object where the data is
field_name : str
name of the field to filter
fix_ranges, fix_azimuths: list of floats
List of ranges [m], azimuths [deg] couples
rng_tol : float
Tolerance between the nominal range and the radar range [m]
ang_tol : float
Tolerance between the nominal angle and the radar angle [deg]
ele_min, ele_max: float
Min and max elevation angle [deg]
Returns
-------
xvals : 1D float array
The ranges of each rng, ele pair
yvals : 1D float array
The values
valid_rng, valid_ele : float arrays
The rng, ele pairs
"""
if ele_min is None:
ele_min = np.min(radar.elevation["data"])
if ele_max is None:
ele_max = np.max(radar.elevation["data"])
yvals = []
xvals = []
valid_rng = []
valid_azi = []
for rng, azi in zip(fix_ranges, fix_azimuths):
ind_rng = find_rng_index(radar.range["data"], rng, rng_tol=rng_tol)
if ind_rng is None:
warn("No range gate found for fix_range " + str(rng))
continue
if radar.scan_type == "ppi":
try:
new_dataset = cross_section_ppi(radar, [azi], az_tol=ang_tol)
except OSError:
warn(
" No data found at range " + str(rng) + " and elevation " + str(azi)
)
continue
else:
ind_sweep = find_ang_index(radar.fixed_angle["data"], azi, ang_tol=ang_tol)
if ind_sweep is None:
warn("No azimuth angle found for fix_azimuth " + str(azi))
continue
new_dataset = radar.extract_sweeps([ind_sweep])
ele_mask = np.logical_and(
new_dataset.elevation["data"] >= ele_min,
new_dataset.elevation["data"] <= ele_max,
)
yvals.extend(new_dataset.fields[field_name]["data"][ele_mask, ind_rng])
xvals.extend(new_dataset.elevation["data"][ele_mask])
valid_rng.append(new_dataset.range["data"][ind_rng])
valid_azi.append(new_dataset.elevation["data"][0])
return np.array(xvals), np.ma.array(yvals), valid_rng, valid_azi
def find_ang_index(ang_vec, ang, ang_tol=0.0):
"""
Find the angle index corresponding to a particular fixed angle
Parameters
----------
ang_vec : float array
The angle data array where to look for
ang : float
The angle to search
ang_tol : float
Tolerance [deg]
Returns
-------
ind_ang : int
The angle index
"""
dist = np.abs(ang_vec - ang)
ind_ang = np.argmin(dist)
if dist[ind_ang] > ang_tol:
return None
return ind_ang
def _create_qvp_object(
radar, field_names, qvp_type="qvp", start_time=None, hmax=10000.0, hres=200.0
):
"""
Creates a QVP object containing fields from a radar object that can
be used to plot and produce the quasi vertical profile
Parameters
----------
radar : Radar
Radar object used.
field_names : list of strings
Radar fields to use for QVP calculation.
qvp_type : str
Type of QVP. Can be qvp, rqvp, evp
start_time : datetime object
the QVP start time
hmax : float
The maximum height of the QVP [m]. Default 10000.
hres : float
The QVP range resolution [m]. Default 50
Returns
-------
qvp : Radar-like object
A quasi vertical profile object containing fields
from a radar object
"""
qvp = deepcopy(radar)
# prepare space for field
qvp.fields = dict()
for field_name in field_names:
qvp.add_field(field_name, deepcopy(radar.fields[field_name]))
qvp.fields[field_name]["data"] = np.array([], dtype="float64")
# fixed radar objects parameters
qvp.range["data"] = np.arange(hmax / hres) * hres + hres / 2.0
qvp.ngates = len(qvp.range["data"])
if start_time is None:
qvp.time["units"] = radar.time["units"]
else:
qvp.time["units"] = make_time_unit_str(start_time)
qvp.time["data"] = np.array([], dtype="float64")
qvp.scan_type = qvp_type
qvp.sweep_number["data"] = np.array([0], dtype="int32")
qvp.nsweeps = 1
qvp.sweep_mode["data"] = np.array([qvp_type])
qvp.sweep_start_ray_index["data"] = np.array([0], dtype="int32")
if qvp.rays_are_indexed is not None:
qvp.rays_are_indexed["data"] = np.array([qvp.rays_are_indexed["data"][0]])
if qvp.ray_angle_res is not None:
qvp.ray_angle_res["data"] = np.array([qvp.ray_angle_res["data"][0]])
if qvp_type in ("rqvp", "evp", "vp"):
qvp.fixed_angle["data"] = np.array([90.0], dtype="float64")
# ray dependent radar objects parameters
qvp.sweep_end_ray_index["data"] = np.array([-1], dtype="int32")
qvp.rays_per_sweep["data"] = np.array([0], dtype="int32")
qvp.azimuth["data"] = np.array([], dtype="float64")
qvp.elevation["data"] = np.array([], dtype="float64")
qvp.nrays = 0
return qvp
def _create_along_coord_object(
radar, field_names, rng_values, fixed_angle, mode, start_time=None
):
"""
Creates an along coord object containing fields from a radar object that
can be used to plot and produce the time series along a coordinate
Parameters
----------
radar : Radar
Radar object used.
field_names : list of strings
Radar fields to use for QVP calculation.
rng_values : 1D-array
The values to put in the range field
fixed_angle : float
the fixed angle
mode : str
The along coord mode, can be ALONG_AZI, ALONG_ELE, ALONG_RNG
start_time : datetime object
the acoord start time
Returns
-------
acoord : Radar-like object
An along coordinate object containing fields from a radar object
"""
acoord = deepcopy(radar)
# prepare space for field
acoord.fields = dict()
for field_name in field_names:
acoord.add_field(field_name, deepcopy(radar.fields[field_name]))
acoord.fields[field_name]["data"] = np.array([], dtype="float64")
# fixed radar objects parameters
acoord.range["data"] = rng_values
acoord.ngates = len(acoord.range["data"])
if start_time is None:
acoord.time["units"] = radar.time["units"]
else:
acoord.time["units"] = make_time_unit_str(start_time)
acoord.time["data"] = np.array([], dtype="float64")
acoord.scan_type = mode
acoord.sweep_number["data"] = np.array([0], dtype="int32")
acoord.nsweeps = 1
acoord.sweep_mode["data"] = np.array([mode])
acoord.sweep_start_ray_index["data"] = np.array([0], dtype="int32")
if acoord.rays_are_indexed is not None:
acoord.rays_are_indexed["data"] = np.array([acoord.rays_are_indexed["data"][0]])
if acoord.ray_angle_res is not None:
acoord.ray_angle_res["data"] = np.array([acoord.ray_angle_res["data"][0]])
acoord.fixed_angle["data"] = np.array([fixed_angle], dtype="float64")
# ray dependent radar objects parameters
acoord.sweep_end_ray_index["data"] = np.array([-1], dtype="int32")
acoord.rays_per_sweep["data"] = np.array([0], dtype="int32")
acoord.azimuth["data"] = np.array([], dtype="float64")
acoord.elevation["data"] = np.array([], dtype="float64")
acoord.nrays = 0
return acoord
def _update_qvp_metadata(qvp, ref_time, lon, lat, elev=90.0):
"""
updates a QVP object metadata with data from the current radar volume
Parameters
----------
qvp : QVP object
QVP object
ref_time : datetime object
the current radar volume reference time
Returns
-------
qvp : QVP object
The updated QVP object
"""
start_time = num2date(0, qvp.time["units"], qvp.time["calendar"])
qvp.time["data"] = np.append(
qvp.time["data"], (ref_time - start_time).total_seconds()
)
qvp.sweep_end_ray_index["data"][0] += 1
qvp.rays_per_sweep["data"][0] += 1
qvp.nrays += 1
qvp.azimuth["data"] = np.ones((qvp.nrays,), dtype="float64") * 0.0
qvp.elevation["data"] = np.ones((qvp.nrays,), dtype="float64") * elev
qvp.gate_longitude["data"] = np.ones((qvp.nrays, qvp.ngates), dtype="float64") * lon
qvp.gate_latitude["data"] = np.ones((qvp.nrays, qvp.ngates), dtype="float64") * lat
qvp.gate_altitude["data"] = ma_broadcast_to(
qvp.range["data"], (qvp.nrays, qvp.ngates)
)
return qvp
def _update_along_coord_metadata(acoord, ref_time, elevation, azimuth):
"""
updates an along coordinate object metadata with data from the current
radar volume
Parameters
----------
acoord : along coordinate object
along coordinate object
ref_time : datetime object
the current radar volume reference time
elevation, azimuth : 1D-array
the elevation and azimuth value of the data selected
Returns
-------
acoord : along coordinate object
The updated along coordinate object
"""
start_time = num2date(0, acoord.time["units"], acoord.time["calendar"])
acoord.time["data"] = np.append(
acoord.time["data"], (ref_time - start_time).total_seconds()
)
acoord.sweep_end_ray_index["data"][0] += 1
acoord.rays_per_sweep["data"][0] += 1
acoord.nrays += 1
acoord.azimuth["data"] = np.append(acoord.azimuth["data"], azimuth)
acoord.elevation["data"] = np.append(acoord.elevation["data"], elevation)
return acoord