Source code for pyart.retrieve.qvp

"""
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)

    # Printing radar tilt angles and radar elevation
    print(radar.fixed_angle["data"])
    print(radar.elevation["data"][-1])

    # 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 in ("rhi", "elevation_surveillance", "manual_rhi"):
        radar_aux = cross_section_rhi(radar_aux, [angle], el_tol=ang_tol)
    elif radar_aux.scan_type in ("ppi", "azimuth_surveillance", "manual_ppi"):
        radar_aux = radar_aux.extract_sweeps([int(angle)])
    else:
        warn(f"Error: unsupported scan type {radar_aux.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 in ("rhi", "elevation_surveillance", "manual_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 in ("ppi", "azimuth_surveillance", "manual_ppi"):
        radar_ppi = radar_aux
    else:
        warn(f"Error: unsupported scan type {radar_aux.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 in ("rhi", "elevation_surveillance", "manual_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 in ("ppi", "azimuth_surveillance", "manual_ppi"):
        radar_ppi = radar_aux
    else:
        warn(f"Error: unsupported scan type {radar_aux.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 in ("rhi", "elevation_surveillance", "manual_rhi"):
        radar_aux = cross_section_rhi(radar_aux, [angle], el_tol=ang_tol)
    elif radar_aux.scan_type in ("ppi", "azimuth_surveillance", "manual_ppi"):
        radar_aux = radar_aux.extract_sweeps([int(angle)])
    else:
        warn(f"Error: unsupported scan type {radar_aux.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 in ("rhi", "elevation_surveillance", "manual_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 in ("ppi", "azimuth_surveillance", "manual_ppi"):
        radar_ppi = radar_aux
    else:
        warn(f"Error: unsupported scan type {radar_aux.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