From d0cb1ba866b3d8aebb1e9da274da89596e9e4fff Mon Sep 17 00:00:00 2001 From: Daniel Wolfensberger Date: Tue, 27 Aug 2024 18:36:03 +0200 Subject: [PATCH 01/15] ENH: addition of qvp functions from MeteoSwiss Py-ART + entropy computation in hydroclass_semisupervised --- pyart/core/__init__.py | 1 + pyart/core/transforms.py | 27 + pyart/retrieve/__init__.py | 3 +- pyart/retrieve/echo_class.py | 514 +++++++--- pyart/retrieve/qpe.py | 141 ++- pyart/retrieve/qvp.py | 1579 ++++++++++++++++++++++++++++- pyart/util/__init__.py | 2 + pyart/util/circular_stats.py | 33 + pyart/util/radar_utils.py | 29 + tests/core/test_transforms.py | 9 + tests/retrieve/test_echo_class.py | 33 +- tests/retrieve/test_qvp.py | 67 ++ tests/util/test_circular_stats.py | 28 + tests/util/test_radar_utils.py | 9 + 14 files changed, 2205 insertions(+), 270 deletions(-) diff --git a/pyart/core/__init__.py b/pyart/core/__init__.py index 15788cac001..abcf1e0c930 100644 --- a/pyart/core/__init__.py +++ b/pyart/core/__init__.py @@ -13,6 +13,7 @@ from .transforms import cartesian_vectors_to_geographic # noqa from .transforms import geographic_to_cartesian # noqa from .transforms import geographic_to_cartesian_aeqd # noqa +from .transforms import cartesian_to_antenna # noqa from .wind_profile import HorizontalWindProfile # noqa __all__ = [s for s in dir() if not s.startswith("_")] diff --git a/pyart/core/transforms.py b/pyart/core/transforms.py index e9581f5c0eb..a4c02f88429 100644 --- a/pyart/core/transforms.py +++ b/pyart/core/transforms.py @@ -77,6 +77,33 @@ def antenna_to_cartesian(ranges, azimuths, elevations): y = s * np.cos(theta_a) return x, y, z +def cartesian_to_antenna(x, y, z): + """ + Returns antenna coordinates from Cartesian coordinates. + + Parameters + ---------- + x, y, z : array + Cartesian coordinates in meters from the radar. + + Returns + ------- + ranges : array + Distances to the center of the radar gates (bins) in m. + azimuths : array + Azimuth angle of the radar in degrees. [-180., 180] + elevations : array + Elevation angle of the radar in degrees. + ke : float, optional + Effective radius scale factor + + """ + ranges = np.sqrt(x ** 2. + y ** 2. + z ** 2.) + elevations = np.rad2deg(np.arctan(z / np.sqrt(x ** 2. + y ** 2.))) + azimuths = np.rad2deg(np.arctan2(x, y)) # [-180, 180] + azimuths[azimuths < 0.] += 360. # [0, 360] + + return ranges, azimuths, elevations def antenna_vectors_to_cartesian(ranges, azimuths, elevations, edges=False): """ diff --git a/pyart/retrieve/__init__.py b/pyart/retrieve/__init__.py index 9be0b022e98..60893b01f29 100644 --- a/pyart/retrieve/__init__.py +++ b/pyart/retrieve/__init__.py @@ -21,7 +21,8 @@ from .qpe import est_rain_rate_zkdp # noqa from .qpe import est_rain_rate_zpoly # noqa from .qpe import ZtoR # noqa -from .qvp import quasi_vertical_profile # noqa +from .qvp import quasi_vertical_profile, compute_qvp, compute_rqvp #noqa +from .qvp import compute_evp, compute_svp, compute_vp, compute_ts_along_coord #noqa from .simple_moment_calculations import calculate_snr_from_reflectivity # noqa from .simple_moment_calculations import calculate_velocity_texture # noqa from .simple_moment_calculations import compute_cdr # noqa diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py index ac19ce5e826..50a3a6d723e 100644 --- a/pyart/retrieve/echo_class.py +++ b/pyart/retrieve/echo_class.py @@ -3,6 +3,7 @@ """ +from copy import deepcopy from warnings import warn import numpy as np @@ -10,6 +11,7 @@ # Local imports from ..config import get_field_name, get_fillvalue, get_metadata from ..core import Grid +from ..util.radar_utils import ma_broadcast_to from ._echo_class import _feature_detection, steiner_class_buff from ._echo_class_wt import calc_scale_break, wavelet_reclass @@ -599,18 +601,19 @@ def feature_detection( return feature_dict -def hydroclass_semisupervised( - radar, - mass_centers=None, - weights=np.array([1.0, 1.0, 1.0, 0.75, 0.5]), - refl_field=None, - zdr_field=None, - rhv_field=None, - kdp_field=None, - temp_field=None, - hydro_field=None, - radar_freq=None, -): +def hydroclass_semisupervised(radar, + hydro_names=('AG', 'CR', 'LR', 'RP', 'RN', 'VI', + 'WS', 'MH', 'IH/HDG'), + var_names=('Zh', 'ZDR', 'KDP', 'RhoHV', + 'relH'), + mass_centers=None, + weights=np.array([1., 1., 1., 0.75, 0.5]), + value=50., lapse_rate=-6.5, refl_field=None, + zdr_field=None, rhv_field=None, kdp_field=None, + temp_field=None, iso0_field=None, + hydro_field=None, entropy_field=None, + temp_ref='temperature', compute_entropy=False, + output_distances=False, vectorize=False): """ Classifies precipitation echoes following the approach by Besic et al (2016). @@ -619,30 +622,47 @@ def hydroclass_semisupervised( ---------- radar : radar Radar object. + hydro_names : array of str + name of the types of hydrometeors + var_names : array of str + name of the variables mass_centers : ndarray 2D, optional The centroids for each variable and hydrometeor class in (nclasses, nvariables). weights : ndarray 1D, optional - The weight given to each variable. - refl_field, zdr_field, rhv_field, kdp_field, temp_field : str, optional + The weight given to each variable. Ordered by [dBZ, ZDR, KDP, RhoHV, + H_ISO0] + value : float + The value controlling the rate of decay in the distance transformation + lapse_rate : float + The decrease in temperature for each vertical km [deg/km] + refl_field, zdr_field, rhv_field, kdp_field, temp_field, iso0_field : str Inputs. Field names within the radar object which represent the horizonal reflectivity, the differential reflectivity, the copolar - correlation coefficient, the specific differential phase and the - temperature field. A value of None for any of these parameters will - use the default field name as defined in the Py-ART configuration - file. - hydro_field : str, optional + correlation coefficient, the specific differential phase, the + temperature (in deg celsius) and the height respect to the iso0 fields. + A value of None for any of these parameters will use the default field + name as defined in the Py-ART configuration file. + hydro_field : str Output. Field name which represents the hydrometeor class field. A value of None will use the default field name as defined in the Py-ART configuration file. - radar_freq : str, optional - Radar frequency in Hertz (Hz) used for classification. - This parameter will be ignored, if the radar object has frequency information. + temp_ref : str + the field use as reference for temperature. Can be either temperature + or height_over_iso0 + compute_entropy : bool + If true, the entropy is computed + output_distances : bool + If true, the normalized distances to the centroids for each + hydrometeor are provided as output + vectorize : bool + If true, a vectorized version of the class assignation is going to be + used Returns ------- - hydro : dict - Hydrometeor classification field. + fields_dict : dict + Dictionary containing the retrieved fields References ---------- @@ -651,99 +671,124 @@ def hydroclass_semisupervised( of polarimetric radar measurements: a semi-supervised approach, Atmos. Meas. Tech., 9, 4425-4445, doi:10.5194/amt-9-4425-2016, 2016 - Notes - ----- - The default hydrometeor classification is valid for C-band radars. For X-band radars, - if frequency information is not present in the `radar.instrument_parameters`, the user-supplied - `radar_freq` will be used with a warning. If both `radar.instrument_parameters` and - `radar_freq` parameter are missing, the algorithm defaults to the C band. - - If the radar frequency information is missing from the radar object, you can add it in - `radar.instrument_parameters`, as follows: - .. code-block:: python - - radar.instrument_parameters["frequency"] = { - "long_name": "Radar frequency", - "units": "Hz", - "data": [9.2e9] - } """ - lapse_rate = -6.5 - # select the centroids as a function of frequency band if mass_centers is None: # assign coefficients according to radar frequency - if radar.instrument_parameters and "frequency" in radar.instrument_parameters: - frequency = radar.instrument_parameters["frequency"]["data"][0] - mass_centers = _get_mass_centers(frequency) - warn(f"Using radar frequency from instrument parameters: {frequency}") - elif radar_freq is not None: - mass_centers = _get_mass_centers(radar_freq) - warn( - f"Radar instrument parameters are empty. Using user-supplied radar frequency: {radar_freq}" - ) + if 'frequency' in radar.instrument_parameters: + mass_centers = _get_mass_centers( + radar.instrument_parameters['frequency']['data'][0]) else: - mass_centers = _mass_centers_table()["C"] - warn( - "Radar instrument parameters and radar_freq param are empty." - "So frequency is unknown. Default coefficients for C band will be applied." - ) + mass_centers = _mass_centers_table()['C'] + warn('Radar frequency unknown. ' + + 'Default coefficients for C band will be applied') - # parse the field parameters - if refl_field is None: - refl_field = get_field_name("reflectivity") - if zdr_field is None: - zdr_field = get_field_name("differential_reflectivity") - if rhv_field is None: - rhv_field = get_field_name("cross_correlation_ratio") - if kdp_field is None: - kdp_field = get_field_name("specific_differential_phase") - if temp_field is None: - temp_field = get_field_name("temperature") if hydro_field is None: - hydro_field = get_field_name("radar_echo_classification") - - # extract fields and parameters from radar - radar.check_field_exists(refl_field) - radar.check_field_exists(zdr_field) - radar.check_field_exists(rhv_field) - radar.check_field_exists(kdp_field) - radar.check_field_exists(temp_field) - - refl = radar.fields[refl_field]["data"] - zdr = radar.fields[zdr_field]["data"] - rhohv = radar.fields[rhv_field]["data"] - kdp = radar.fields[kdp_field]["data"] - temp = radar.fields[temp_field]["data"] - - # convert temp in relative height respect to iso0 - relh = temp * (1000.0 / lapse_rate) - - # standardize data - refl_std = _standardize(refl, "Zh") - zdr_std = _standardize(zdr, "ZDR") - kdp_std = _standardize(kdp, "KDP") - rhohv_std = _standardize(rhohv, "RhoHV") - relh_std = _standardize(relh, "relH") - - # standardize centroids - mc_std = np.zeros(np.shape(mass_centers)) - mc_std[:, 0] = _standardize(mass_centers[:, 0], "Zh") - mc_std[:, 1] = _standardize(mass_centers[:, 1], "ZDR") - mc_std[:, 2] = _standardize(mass_centers[:, 2], "KDP") - mc_std[:, 3] = _standardize(mass_centers[:, 3], "RhoHV") - mc_std[:, 4] = _standardize(mass_centers[:, 4], "relH") + hydro_field = get_field_name('radar_echo_classification') + if compute_entropy: + if entropy_field is None: + entropy_field = get_field_name('hydroclass_entropy') + + # Get the data fields + fields_dict = {} + for var_name in var_names: + if var_name == 'Zh': + if refl_field is None: + refl_field = get_field_name('reflectivity') + radar.check_field_exists(refl_field) + fields_dict.update({ + var_name: deepcopy(radar.fields[refl_field]['data'])}) + elif var_name == 'ZDR': + if zdr_field is None: + zdr_field = get_field_name('differential_reflectivity') + radar.check_field_exists(zdr_field) + fields_dict.update({ + var_name: deepcopy(radar.fields[zdr_field]['data'])}) + elif var_name == 'KDP': + if kdp_field is None: + kdp_field = get_field_name('specific_differential_phase') + radar.check_field_exists(kdp_field) + fields_dict.update({ + var_name: deepcopy(radar.fields[kdp_field]['data'])}) + elif var_name == 'RhoHV': + if rhv_field is None: + rhv_field = get_field_name('cross_correlation_ratio') + radar.check_field_exists(rhv_field) + fields_dict.update({ + var_name: deepcopy(radar.fields[rhv_field]['data'])}) + elif var_name == 'relH': + if temp_ref == 'temperature': + if temp_field is None: + temp_field = get_field_name('temperature') + radar.check_field_exists(temp_field) + # convert temp in relative height respect to iso0 + temp = deepcopy(radar.fields[temp_field]['data']) + fields_dict.update({var_name: temp * (1000. / lapse_rate)}) + else: + if iso0_field is None: + iso0_field = get_field_name('height_over_iso0') + radar.check_field_exists(iso0_field) + fields_dict.update({ + var_name: deepcopy(radar.fields[iso0_field]['data'])}) + else: + raise ValueError( + 'Variable ' + var_name + ' unknown. ' + + 'Valid variable names for hydrometeor classification are: ' + + 'relH, Zh, ZDR, KDP and RhoHV') + + # standardize data and centroids + mc_std = np.empty( + np.shape(mass_centers), dtype=fields_dict[var_names[0]].dtype) + for i, var_name in enumerate(var_names): + mc_std[:, i] = _standardize(mass_centers[:, i], var_name) + fields_dict[var_name] = _standardize(fields_dict[var_name], var_name) + + # if entropy has to be computed get transformation parameters + t_vals = None + if compute_entropy: + t_vals = _compute_coeff_transform( + mc_std, weights=weights, value=value) # assign to class - hydroclass_data, min_dist = _assign_to_class( - refl_std, zdr_std, kdp_std, rhohv_std, relh_std, mc_std, weights=weights - ) + if vectorize: + hydroclass_data, entropy_data, prop_data = _assign_to_class_scan( + fields_dict, mc_std, var_names=var_names, weights=weights, + t_vals=t_vals) + else: + hydroclass_data, entropy_data, prop_data = _assign_to_class( + fields_dict, mc_std, weights=weights, t_vals=t_vals) # prepare output fields + fields_dict = dict() hydro = get_metadata(hydro_field) - hydro["data"] = hydroclass_data - - return hydro + hydro['data'] = hydroclass_data + hydro.update({'_FillValue': 0}) + labels = ['NC'] + ticks = [1] + boundaries = [0.5, 1.5] + for i, hydro_name in enumerate(hydro_names): + labels.append(hydro_name) + ticks.append(i + 2) + boundaries.append(i + 2.5) + hydro.update({ + 'labels': labels, + 'ticks': ticks, + 'boundaries': boundaries}) + fields_dict.update({'hydro': hydro}) + + if compute_entropy: + entropy = get_metadata(entropy_field) + entropy['data'] = entropy_data + fields_dict.update({'entropy': entropy}) + + if output_distances: + for field_name in hydro_names: + field_name = 'proportion_' + field_name + prop = get_metadata(field_name) + prop['data'] = prop_data[:, :, i] + fields_dict.update({field_name: prop}) + + return fields_dict def _standardize(data, field_name, mx=None, mn=None): @@ -795,74 +840,203 @@ def _standardize(data, field_name, mx=None, mn=None): return field_std - -def _assign_to_class( - zh, - zdr, - kdp, - rhohv, - relh, - mass_centers, - weights=np.array([1.0, 1.0, 1.0, 0.75, 0.5]), -): +def _assign_to_class(fields_dict, mass_centers, + var_names=('Zh', 'ZDR', 'KDP', 'RhoHV', 'relH'), + weights=np.array([1., 1., 1., 0.75, 0.5]), + t_vals=None): """ Assigns an hydrometeor class to a radar range bin computing the distance between the radar variables an a centroid. Parameters ---------- - zh, zdr, kdp, rhohv, relh : radar fields - Variables used for assignment normalized to [-1, 1] values. + fields_dict : dict + Dictionary containg the variables used for assigment normalized to + [-1, 1] values + mass_centers : matrix + centroids normalized to [-1, 1] values (nclasses, nvariables) + var_names : array of str + Name of the variables + weights : array + optional. The weight given to each variable (nvariables) + t_vals : array + transformation values for the distance to centroids (nclasses) + + Returns + ------- + hydroclass : int array + the index corresponding to the assigned class + entropy : float array + the entropy + t_dist : float matrix + if entropy is computed, the transformed distances of each class + (proxy for proportions of each hydrometeor) (nrays, nbins, nclasses) + + """ + # prepare data + nrays = fields_dict[var_names[0]].shape[0] + nbins = fields_dict[var_names[0]].shape[1] + nclasses = mass_centers.shape[0] + nvariables = mass_centers.shape[1] + dtype = fields_dict[var_names[0]].dtype + + hydroclass = np.ma.empty((nrays, nbins), dtype=np.uint8) + entropy = None + t_dist = None + if t_vals is not None: + entropy = np.ma.empty((nrays, nbins), dtype=dtype) + t_dist = np.ma.masked_all((nrays, nbins, nclasses), dtype=dtype) + + for ray in range(nrays): + data = [] + for var_name in var_names: + data.append(fields_dict[var_name][ray, :]) + data = np.ma.array(data, dtype=dtype) + weights_mat = np.broadcast_to( + weights.reshape(nvariables, 1), (nvariables, nbins)) + dist = np.ma.zeros((nclasses, nbins), dtype=dtype) + + # compute distance: masked entries will not contribute to the distance + mask = np.ma.getmaskarray(fields_dict[var_names[0]][ray, :]) + for i in range(nclasses): + centroids_class = mass_centers[i, :] + centroids_class = np.broadcast_to( + centroids_class.reshape(nvariables, 1), (nvariables, nbins)) + dist_ray = np.ma.sqrt(np.ma.sum( + ((centroids_class - data)**2.) * weights_mat, axis=0)) + dist_ray[mask] = np.ma.masked + dist[i, :] = dist_ray + + # Get hydrometeor class + class_vec = dist.argsort(axis=0, fill_value=10e40) + hydroclass_ray = (class_vec[0, :] + 2).astype(np.uint8) + hydroclass_ray[mask] = 1 + hydroclass[ray, :] = hydroclass_ray + + if t_vals is None: + continue + + # Transform the distance using the coefficient of the dominant class + t_vals_ray = np.ma.masked_where(mask, t_vals[class_vec[0, :]]) + t_vals_ray = ma_broadcast_to( + t_vals_ray.reshape(1, nbins), (nclasses, nbins)) + t_dist_ray = np.ma.exp(-t_vals_ray * dist) + + # set transformed distances to a value between 0 and 1 + dist_total = np.ma.sum(t_dist_ray, axis=0) + dist_total = ma_broadcast_to( + dist_total.reshape(1, nbins), (nclasses, nbins)) + t_dist_ray /= dist_total + + # Compute entropy + entropy_ray = -np.ma.sum( + t_dist_ray * np.ma.log(t_dist_ray) / np.ma.log(nclasses), axis=0) + entropy_ray[mask] = np.ma.masked + entropy[ray, :] = entropy_ray + + t_dist[ray, :, :] = np.ma.transpose(t_dist_ray) + + if t_vals is not None: + t_dist *= 100. + + return hydroclass, entropy, t_dist + + +def _assign_to_class_scan(fields_dict, mass_centers, + var_names=('Zh', 'ZDR', 'KDP', 'RhoHV', 'relH'), + weights=np.array([1., 1., 1., 0.75, 0.5]), + t_vals=None): + """ + assigns an hydrometeor class to a radar range bin computing + the distance between the radar variables an a centroid. + Computes the entire radar volume at once + + Parameters + ---------- + fields_dict : dict + Dictionary containg the variables used for assigment normalized to + [-1, 1] values mass_centers : matrix - Centroids normalized to [-1, 1] values. - weights : array, optional - The weight given to each variable. + centroids normalized to [-1, 1] values + var_names : array of str + Name of the variables + weights : array + optional. The weight given to each variable + t_vals : matrix + transformation values for the distance to centroids + (nclasses, nvariables) Returns ------- hydroclass : int array - The index corresponding to the assigned class. - mind_dist : float array - The minimum distance to the centroids. + the index corresponding to the assigned class + entropy : float array + the entropy + t_dist : float matrix + if entropy is computed, the transformed distances of each class + (proxy for proportions of each hydrometeor) (nrays, nbins, nclasses) """ # prepare data - nrays = zh.shape[0] - nbins = zdr.shape[1] + nrays = fields_dict[var_names[0]].shape[0] + nbins = fields_dict[var_names[0]].shape[1] nclasses = mass_centers.shape[0] nvariables = mass_centers.shape[1] + dtype = fields_dict[var_names[0]].dtype - data = np.ma.array([zh, zdr, kdp, rhohv, relh]) + data = [] + for var_name in var_names: + data.append(fields_dict[var_name]) + data = np.ma.array(data, dtype=dtype) weights_mat = np.broadcast_to( - weights.reshape(nvariables, 1, 1), (nvariables, nrays, nbins) - ) - dist = np.ma.zeros((nclasses, nrays, nbins), dtype="float64") + weights.reshape(nvariables, 1, 1), (nvariables, nrays, nbins)) # compute distance: masked entries will not contribute to the distance + mask = np.ma.getmaskarray(fields_dict[var_names[0]]) + dist = np.ma.zeros((nrays, nbins, nclasses), dtype=dtype) + t_dist = None + entropy = None for i in range(nclasses): centroids_class = mass_centers[i, :] centroids_class = np.broadcast_to( - centroids_class.reshape(nvariables, 1, 1), (nvariables, nrays, nbins) - ) - dist[i, :, :] = np.ma.sqrt( - np.ma.sum(((centroids_class - data) ** 2.0) * weights_mat, axis=0) - ) - - # use very large fill_value so that masked entries will be sorted at the - # end. There should not be any masked entry anyway - class_vec = dist.argsort(axis=0, fill_value=10e40) - - # get minimum distance. Acts as a confidence value - dist.sort(axis=0, fill_value=10e40) - min_dist = dist[0, :, :] - - # Entries with non-valid reflectivity values are set to 0 (No class) - mask = np.ma.getmaskarray(zh) - hydroclass = class_vec[0, :, :] + 1 - hydroclass[mask] = 0 - - return hydroclass, min_dist - + centroids_class.reshape(nvariables, 1, 1), + (nvariables, nrays, nbins)) + dist_aux = np.ma.sqrt(np.ma.sum( + ((centroids_class - data)**2.) * weights_mat, axis=0)) + dist_aux[mask] = np.ma.masked + dist[:, :, i] = dist_aux + + del data + del weights_mat + + # Get hydrometeor class + class_vec = dist.argsort(axis=-1, fill_value=10e40) + hydroclass = np.ma.asarray(class_vec[:, :, 0] + 2, dtype=np.uint8) + hydroclass[mask] = 1 + + if t_vals is not None: + # Transform the distance using the coefficient of the dominant class + t_vals_aux = np.ma.masked_where(mask, t_vals[class_vec[:, :, 0]]) + t_vals_aux = ma_broadcast_to( + t_vals_aux.reshape(nrays, nbins, 1), (nrays, nbins, nclasses)) + t_dist = np.ma.exp(-t_vals_aux * dist) + del t_vals_aux + + # set distance to a value between 0 and 1 + dist_total = np.ma.sum(t_dist, axis=-1) + dist_total = ma_broadcast_to( + dist_total.reshape(nrays, nbins, 1), (nrays, nbins, nclasses)) + t_dist /= dist_total + del dist_total + + # compute entroy + entropy = -np.ma.sum( + t_dist * np.ma.log(t_dist) / np.ma.log(nclasses), axis=-1) + entropy[mask] = np.ma.masked + + t_dist *= 100. + + return hydroclass, entropy, t_dist def _get_mass_centers(freq): """ @@ -1165,3 +1339,43 @@ def conv_strat_raut( } return reclass_dict + +def _compute_coeff_transform(mass_centers, + weights=np.array([1., 1., 1., 0.75, 0.5]), + value=50.): + """ + get the transformation coefficients + + Parameters + ---------- + mass_centers : ndarray 2D + The centroids for each class and variable (nclasses, nvariables) + weights : array + optional. The weight given to each variable (nvariables) + value : float + parameter controlling the rate of decay of the distance transformation + + Returns + ------- + t_vals : ndarray 1D + The coefficients used to transform the distances to each centroid for + each class (nclasses) + + """ + nclasses, nvariables = np.shape(mass_centers) + t_vals = np.empty((nclasses, nclasses), dtype=mass_centers.dtype) + for i in range(nclasses): + weights_mat = np.broadcast_to( + weights.reshape(1, nvariables), (nclasses, nvariables)) + centroids_class = mass_centers[i, :] + centroids_class = np.broadcast_to( + centroids_class.reshape(1, nvariables), (nclasses, nvariables)) + t_vals[i, :] = np.sqrt( + np.sum(weights_mat * np.power( + np.abs(centroids_class - mass_centers), 2.), axis=1)) + + # pick the second lowest value (the first is 0) + t_vals = np.sort(t_vals, axis=-1)[:, 1] + t_vals = np.log(value) / t_vals + + return t_vals diff --git a/pyart/retrieve/qpe.py b/pyart/retrieve/qpe.py index c657228d826..27bfd9ca5f2 100644 --- a/pyart/retrieve/qpe.py +++ b/pyart/retrieve/qpe.py @@ -416,24 +416,11 @@ def est_rain_rate_za( return rain_main - -def est_rain_rate_hydro( - radar, - alphazr=0.0376, - betazr=0.6112, - alphazs=0.1, - betazs=0.5, - alphaa=None, - betaa=None, - mp_factor=0.6, - refl_field=None, - a_field=None, - hydro_field=None, - rr_field=None, - main_field=None, - thresh=None, - thresh_max=False, -): +def est_rain_rate_hydro(radar, alphazr=0.0376, betazr=0.6112, alphazs=0.1, + betazs=0.5, alphaa=None, betaa=None, mp_factor=0.6, + refl_field=None, a_field=None, hydro_field=None, + rr_field=None, master_field=None, thresh=None, + thresh_max=False): """ Estimates rainfall rate using different relations between R and the polarimetric variables depending on the hydrometeor type. @@ -460,15 +447,15 @@ def est_rain_rate_hydro( Name of the hydrometeor classification field to use. rr_field : str, optional Name of the rainfall rate field. - main_field : str, optional - Name of the field that is going to act as main. Has to be + master_field : str, optional + Name of the field that is going to act as master. Has to be either refl_field or kdp_field. Default is refl_field. thresh : float, optional - Value of the threshold that determines when to use the secondary + Value of the threshold that determines when to use the slave field. thresh_max : Bool, optional - If true the main field is used up to the thresh value maximum. - Otherwise the main field is not used below thresh value. + If true the master field is used up to the thresh value maximum. + Otherwise the master field is not used below thresh value. Returns ------- @@ -478,90 +465,98 @@ def est_rain_rate_hydro( """ # parse the field parameters if refl_field is None: - refl_field = get_field_name("reflectivity") + refl_field = get_field_name('reflectivity') if a_field is None: - a_field = get_field_name("specific_attenuation") + a_field = get_field_name('specific_attenuation') if hydro_field is None: - hydro_field = get_field_name("radar_echo_classification") + hydro_field = get_field_name('radar_echo_classification') if rr_field is None: - rr_field = get_field_name("radar_estimated_rain_rate") + rr_field = get_field_name('radar_estimated_rain_rate') + if master_field is None: + master_field = refl_field # extract fields and parameters from radar if hydro_field in radar.fields: - hydroclass = radar.fields[hydro_field]["data"] + hydroclass = radar.fields[hydro_field]['data'] else: - raise KeyError("Field not available: " + hydro_field) + raise KeyError('Field not available: ' + hydro_field) # get the location of each hydrometeor class - is_ds = hydroclass == 1 - is_cr = hydroclass == 2 - is_lr = hydroclass == 3 - is_gr = hydroclass == 4 - is_rn = hydroclass == 5 - is_vi = hydroclass == 6 - is_ws = hydroclass == 7 - is_mh = hydroclass == 8 - is_ih = hydroclass == 9 + is_ds = hydroclass == 2 + is_cr = hydroclass == 3 + is_lr = hydroclass == 4 + is_gr = hydroclass == 5 + is_rn = hydroclass == 6 + is_vi = hydroclass == 7 + is_ws = hydroclass == 8 + is_mh = hydroclass == 9 + is_ih = hydroclass == 10 # compute z-r (in rain) z-r in snow and z-a relations rain_z = est_rain_rate_z( - radar, alpha=alphazr, beta=betazr, refl_field=refl_field, rr_field=rr_field - ) + radar, alpha=alphazr, beta=betazr, refl_field=refl_field, + rr_field=rr_field) snow_z = est_rain_rate_z( - radar, alpha=alphazs, beta=betazs, refl_field=refl_field, rr_field=rr_field - ) + radar, alpha=alphazs, beta=betazs, refl_field=refl_field, + rr_field=rr_field) rain_a = est_rain_rate_a( - radar, alpha=alphaa, beta=betaa, a_field=a_field, rr_field=rr_field - ) + radar, alpha=alphaa, beta=betaa, a_field=a_field, rr_field=rr_field) # initialize rainfall rate field - rr_data = np.ma.zeros(hydroclass.shape, dtype="float32") + rr_data = np.ma.zeros(hydroclass.shape, dtype='float32') rr_data[:] = np.ma.masked rr_data.set_fill_value(get_fillvalue()) # apply the relations for each type # solid phase - rr_data[is_ds] = snow_z["data"][is_ds] - rr_data[is_cr] = snow_z["data"][is_cr] - rr_data[is_vi] = snow_z["data"][is_vi] - rr_data[is_gr] = snow_z["data"][is_gr] - rr_data[is_ih] = snow_z["data"][is_ih] + rr_data[is_ds] = snow_z['data'][is_ds] + rr_data[is_cr] = snow_z['data'][is_cr] + rr_data[is_vi] = snow_z['data'][is_vi] + rr_data[is_gr] = snow_z['data'][is_gr] + rr_data[is_ih] = snow_z['data'][is_ih] # rain - if main_field == refl_field: - rain_main = rain_z - rain_secondary = rain_a - elif main_field == a_field: - rain_main = rain_a - rain_secondary = rain_z - elif main_field is None: - main_field = a_field - rain_main = rain_a - rain_secondary = rain_z + if master_field == refl_field: + rain_master = rain_z + rain_slave = rain_a + if thresh is None: + thresh = 5. + thresh_max = True + elif master_field == a_field: + rain_master = rain_a + rain_slave = rain_z + if thresh is None: + thresh = 5. + thresh_max = False + elif master_field is None: + rain_master = rain_a + rain_slave = rain_z + thresh = 5. + thresh_max = False else: - main_field = a_field - rain_main = rain_a - rain_secondary = rain_z - thresh = 0.04 + rain_master = rain_a + rain_slave = rain_z + thresh = 5. thresh_max = False - warn("Unknown main field. Using " + a_field + " with threshold " + str(thresh)) + warn('Unknown master field. Using ' + a_field + ' with threshold ' + + str(thresh)) if thresh_max: - is_secondary = rain_main["data"] > thresh + is_slave = rain_master['data'] > thresh else: - is_secondary = rain_main["data"] < thresh + is_slave = rain_master['data'] < thresh - rain_main["data"][is_secondary] = rain_secondary["data"][is_secondary] + rain_master['data'][is_slave] = rain_slave['data'][is_slave] - rr_data[is_lr] = rain_main["data"][is_lr] - rr_data[is_rn] = rain_main["data"][is_rn] + rr_data[is_lr] = rain_master['data'][is_lr] + rr_data[is_rn] = rain_master['data'][is_rn] # mixed phase - rr_data[is_ws] = mp_factor * rain_z["data"][is_ws] - rr_data[is_mh] = mp_factor * rain_z["data"][is_mh] + rr_data[is_ws] = mp_factor * rain_z['data'][is_ws] + rr_data[is_mh] = mp_factor * rain_z['data'][is_mh] rain = get_metadata(rr_field) - rain["data"] = rr_data + rain['data'] = rr_data return rain diff --git a/pyart/retrieve/qvp.py b/pyart/retrieve/qvp.py index 991b43d5930..562f1236eb6 100644 --- a/pyart/retrieve/qvp.py +++ b/pyart/retrieve/qvp.py @@ -1,16 +1,48 @@ """ +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 -def quasi_vertical_profile( - radar, desired_angle=None, fields=None, gatefilter=None, verbose=False -): +def quasi_vertical_profile(radar, desired_angle=None, fields=None, + gatefilter=None): """ Quasi Vertical Profile. @@ -26,8 +58,6 @@ def quasi_vertical_profile( desired_angle : float Radar tilt angle to use for indexing radar field data. None will result in wanted_angle = 20.0 - verbose : bool - boolean flag to turn the printing of radar tilts on or off. Other Parameters ---------------- @@ -38,38 +68,38 @@ def quasi_vertical_profile( Returns ------- qvp : Dictonary - A quasai vertical profile object containing fields + 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.. + 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. + 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. + 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 + # 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() + index = abs(radar.fixed_angle['data'] - desired_angle).argmin() radar_slice = radar.get_slice(index) - if verbose: - # Printing radar tilt angles and radar elevation - print(radar.fixed_angle["data"]) - print(radar.elevation["data"][-1]) + # 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 @@ -77,13 +107,13 @@ def quasi_vertical_profile( 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 - ) + 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) @@ -91,13 +121,12 @@ def quasi_vertical_profile( qvp.update({field: radar_fields}) else: - # Filtereing data based on defined gatefilter + # 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 - ) + 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) @@ -105,9 +134,1503 @@ def quasi_vertical_profile( 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}) + 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 + + +def compute_qvp(radar, field_names, ref_time=None, angle=0., ang_tol=1., + hmax=10000., hres=50., 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 + + +def compute_rqvp(radar, field_names, ref_time=None, hmax=10000., hres=2., + avg_type='mean', nvalid_min=30, interp_kind='nearest', + rmax=50000., weight_power=2., 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.) + + rmax_km = rmax / 1000. + 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.) + + np.power(radar_aux.gate_y['data'][0, :], 2.)) / 1000. + + # 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.)) + weight[grng_interp <= rmax_km - 1.] = 1. / np.power( + weight[grng_interp <= rmax_km - 1.], 0.) + + if weight_power == -1: + weight[grng_interp > rmax_km - 1.] = 0. + else: + weight[grng_interp > rmax_km - 1.] = 1. / np.power( + weight[grng_interp > rmax_km - 1.], 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 + + +def compute_evp(radar, field_names, lon, lat, ref_time=None, + latlon_tol=0.0005, delta_rng=15000., delta_azi=10, + hmax=10000., hres=250., 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.) + + 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 + + +def compute_svp(radar, field_names, lon, lat, angle, ref_time=None, + ang_tol=1., latlon_tol=0.0005, delta_rng=15000., delta_azi=10, + hmax=10000., hres=250., 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 + + +def compute_vp(radar, field_names, lon, lat, ref_time=None, + latlon_tol=0.0005, hmax=10000., hres=50., + 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.) + + 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 + + +def compute_ts_along_coord(radar, field_names, mode='ALONG_AZI', + fixed_range=None, fixed_azimuth=None, + fixed_elevation=None, ang_tol=1., rng_tol=50., + 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. + if value_stop is None: + value_stop = radar.range['data'][-1] + if fixed_azimuth is None: + fixed_azimuth = 0. + if fixed_elevation is None: + fixed_elevation = 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. + if fixed_elevation is None: + fixed_elevation = 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. + if fixed_azimuth is None: + fixed_azimuth = 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.): + """ + 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.) + 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.): + """ + 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.: + azi_max -= 360. + if azi_min < 0.: + azi_min += 360. + 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., 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. + 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., ang_tol=1., 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., ang_tol=1., 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.): + """ + 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., hres=200.): + """ + 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. + 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.], 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.): + """ + 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. + 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 diff --git a/pyart/util/__init__.py b/pyart/util/__init__.py index ee8e1ae47cd..1a1674534de 100644 --- a/pyart/util/__init__.py +++ b/pyart/util/__init__.py @@ -14,6 +14,7 @@ from .circular_stats import interval_std # noqa from .circular_stats import mean_of_two_angles # noqa from .circular_stats import mean_of_two_angles_deg # noqa +from .circular_stats import compute_directional_stats # noqa from .columnsect import for_azimuth # noqa from .columnsect import get_column_rays # noqa from .columnsect import get_field_location # noqa @@ -32,6 +33,7 @@ determine_sweeps, subset_radar, to_vpt, + ma_broadcast_to ) from .sigmath import ( # noqa angular_texture_2d, diff --git a/pyart/util/circular_stats.py b/pyart/util/circular_stats.py index 639cbb83b26..b9a6003e1cb 100644 --- a/pyart/util/circular_stats.py +++ b/pyart/util/circular_stats.py @@ -9,6 +9,39 @@ # https://en.wikipedia.org/wiki/Directional_statistics # https://en.wikipedia.org/wiki/Mean_of_circular_quantities +def compute_directional_stats(field, avg_type='mean', nvalid_min=1, axis=0): + """ + Computes the mean or the median along one of the axis (ray or range) + + Parameters + ---------- + field : ndarray + the radar field data + avg_type :str + the type of average: 'mean' or 'median' + nvalid_min : int + the minimum number of points to consider the stats valid. Default 1 + axis : int + the axis along which to compute (0=ray, 1=range) + + Returns + ------- + values : ndarray 1D + The resultant statistics + nvalid : ndarray 1D + The number of valid points used in the computation + """ + if avg_type == 'mean': + values = np.ma.mean(field, axis=axis) + else: + values = np.ma.median(field, axis=axis) + + # Set to non-valid if there is not a minimum number of valid gates + valid = np.logical_not(np.ma.getmaskarray(field)) + nvalid = np.sum(valid, axis=axis, dtype=int) + values[nvalid < nvalid_min] = np.ma.masked + + return values, nvalid def mean_of_two_angles(angles1, angles2): """ diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index a8f056728c6..482068d022c 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -689,3 +689,32 @@ def image_mute_radar(radar, field, mute_field, mute_threshold, field_threshold=N muted_dict["long_name"] = "Muted " + field radar.add_field("muted_" + field, muted_dict) return radar + +def ma_broadcast_to(array, tup): + """ + Is used to guarantee that a masked array can be broadcasted without + loosing the mask + + Parameters + ---------- + array : Numpy masked array or normal array + tup : shape as tuple + + Returns + ------- + broadcasted_array + The broadcasted numpy array including its mask if available + otherwise only the broadcasted array is returned + + """ + broadcasted_array = np.broadcast_to(array, tup) + + if np.ma.is_masked(array): + initial_mask = np.ma.getmask(array) + initial_fill_value = array.fill_value + broadcasted_mask = np.broadcast_to(initial_mask, tup) + return np.ma.array( + broadcasted_array, mask=broadcasted_mask, + fill_value=initial_fill_value) + + return broadcasted_array \ No newline at end of file diff --git a/tests/core/test_transforms.py b/tests/core/test_transforms.py index ab7527bbf39..e07124f4337 100644 --- a/tests/core/test_transforms.py +++ b/tests/core/test_transforms.py @@ -188,3 +188,12 @@ def test_cartesian_to_geographic_aeqd(): lon, lat = transforms.cartesian_to_geographic_aeqd(x, y, lon_0, lat_0, R) assert_almost_equal(lon, -100.0, 3) assert_almost_equal(lat, 40.0, 3) + +def test_cartesian_to_antenna(): + r,az,el = transforms.cartesian_to_antenna(np.array([1000,3000,2000]), + np.array([1000,2000,3000]), + np.array([500,1000,1500])) + + assert np.allclose(r, [1500, 3741.6573, 3905.1248]) + assert np.allclose(az, [45., 56.30993247, 33.69006753]) + assert np.allclose(el, [19.47122063, 15.50135957, 22.5885387]) diff --git a/tests/retrieve/test_echo_class.py b/tests/retrieve/test_echo_class.py index ce9e276250e..fbe3ecf6f36 100644 --- a/tests/retrieve/test_echo_class.py +++ b/tests/retrieve/test_echo_class.py @@ -144,27 +144,28 @@ def test_hydroclass_semisupervised(): mass_centers = pyart.retrieve.echo_class._get_mass_centers(10e9) - hydro_nofreq = pyart.retrieve.hydroclass_semisupervised(radar_small) + hydro_nofreq = pyart.retrieve.hydroclass_semisupervised(radar_small)['hydro'] + assert "units" in hydro_nofreq.keys() assert "standard_name" in hydro_nofreq.keys() assert "long_name" in hydro_nofreq.keys() assert "coordinates" in hydro_nofreq.keys() - assert hydro_nofreq["data"].dtype == "int64" - assert_allclose(hydro_nofreq["data"][0][0:5], [6, 6, 6, 6, 6]) - assert_allclose(hydro_nofreq["data"][0][-5:], [2, 2, 2, 2, 2]) + assert hydro_nofreq["data"].dtype == "uint8" + assert_allclose(hydro_nofreq["data"][0][0:5], [7, 7, 7, 7, 7]) + assert_allclose(hydro_nofreq["data"][0][-5:], [3, 3, 3, 3, 3]) radar_small.instrument_parameters["frequency"] = {"data": np.array([5e9])} - hydro_freq = pyart.retrieve.hydroclass_semisupervised(radar_small) - assert_allclose(hydro_freq["data"][0][0:5], [6, 6, 6, 6, 6]) - assert_allclose(hydro_freq["data"][0][-5:], [2, 2, 2, 2, 2]) + hydro_freq = pyart.retrieve.hydroclass_semisupervised(radar_small)['hydro'] + assert_allclose(hydro_freq["data"][0][0:5], [7, 7, 7, 7, 7]) + assert_allclose(hydro_freq["data"][0][-5:], [3, 3, 3, 3, 3]) hydro = pyart.retrieve.hydroclass_semisupervised( radar_small, mass_centers=mass_centers - ) - assert_allclose(hydro["data"][0][0:5], [6, 6, 6, 6, 6]) - assert_allclose(hydro["data"][0][-5:], [2, 2, 2, 2, 2]) + )['hydro'] + assert_allclose(hydro["data"][0][0:5], [7, 7, 7, 7, 7]) + assert_allclose(hydro["data"][0][-5:], [3, 3, 3, 3, 3]) def test_data_limits_table(): dlimits_dict = pyart.retrieve.echo_class._data_limits_table() @@ -259,16 +260,12 @@ def test_assign_to_class(): ) mass_centers = pyart.retrieve.echo_class._get_mass_centers(10e9) - hydroclass, min_dist = pyart.retrieve.echo_class._assign_to_class( - zh, zdr, kdp, rhohv, relh, mass_centers - ) + field_dict = {'Zh':zh, 'ZDR':zdr, 'KDP':kdp, 'RhoHV':rhohv, 'relH':relh} - assert_allclose(hydroclass[0], [7, 7, 7, 7, 7], atol=1e-7) - assert_allclose( - min_dist[0], - [227.0343910, 227.0343910, 227.0343910, 227.0343910, 227.0343910], - atol=1e-7, + hydroclass, _, _ = pyart.retrieve.echo_class._assign_to_class( + field_dict, mass_centers ) + assert_allclose(hydroclass[0], [8, 8, 8, 8, 8], atol=1e-7) def test_standardize(): diff --git a/tests/retrieve/test_qvp.py b/tests/retrieve/test_qvp.py index d292dfa4ffd..211d62b0ac3 100644 --- a/tests/retrieve/test_qvp.py +++ b/tests/retrieve/test_qvp.py @@ -1,11 +1,78 @@ """ Unit Tests for Py-ART's retrieve/qvp.py module. """ import numpy as np +import pytest from numpy.testing import assert_almost_equal import pyart +@pytest.fixture +def test_radar(): + test_radar = pyart.testing.make_empty_ppi_radar(1000, 360, 1) + test_radar.range['data'] *= 100 + test_radar.elevation['data'] = np.ones(test_radar.elevation['data'].shape) * 10. + test_radar.fixed_angle['data'] = np.array([10.]) + refl = 0.1 * np.arange(test_radar.ngates) + refl = np.tile(refl, (test_radar.nrays, 1)) + test_radar.add_field('reflectivity', {'data': refl}) + return test_radar + +def test_compute_qvp(test_radar): + qvp = pyart.retrieve.compute_qvp(test_radar, ['reflectivity'], hmax = 10000) + qvp_refl = np.array([1.8999999999999875, 2.200000000000011, + 2.3999999999999853, 2.7000000000000175, 3.0, + 3.299999999999978, 3.599999999999994, + 3.9000000000000132, 4.200000000000027, 4.5, + 4.7000000000000295, 5.0, 5.299999999999969, + 5.599999999999963, 5.900000000000039, + 6.2000000000000135, 6.5, 6.7, 7.0, 7.300000000000015]) + + assert np.allclose(qvp.fields['reflectivity']['data'][0,10:30], qvp_refl) + assert len(qvp.range['data']) == 200 + assert np.all(qvp.range['data'] == np.arange(25, 10000, 50)) + assert qvp.azimuth['data'][0] == 0 + +def test_compute_rqvp(test_radar): + rqvp = pyart.retrieve.compute_rqvp(test_radar, ['reflectivity'], hres = 50, + hmax = 5000, rmax = 5000, weight_power = -1) + rqvp_refl = np.array([1.8999999999999875, 2.200000000000011, + 2.3999999999999853, 2.7000000000000175, 3.0, + 3.299999999999978, 3.599999999999994, + 3.9000000000000132, np.nan, np.nan, np.nan, np.nan, + np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, + np.nan]) + rqvp_refl_mask = np.array([False, False, False, False, False, False, False, False, + True, True, True, True, True, True, True, True, + True, True, True, True]) + + assert np.allclose(rqvp.fields['reflectivity']['data'][0,10:30], rqvp_refl) + assert np.allclose(rqvp.fields['reflectivity']['data'][0,10:30].mask, rqvp_refl_mask) + + assert len(rqvp.range['data']) == 100 + assert np.all(rqvp.range['data'] == np.arange(25, 5000, 50)) + assert rqvp.azimuth['data'][0] == 0 + +def test_compute_evp(test_radar): + evp = pyart.retrieve.compute_evp(test_radar, ['reflectivity'], -97, 36, + delta_rng = 40000) + + evp_refl = np.array([np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, + np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, + 32.7, 33.7, 35.10000000000001, 36.5, 37.89999999999999, + 39.29999999999999, 40.70000000000001]) + + evp_refl_mask = np.array([[ True, True, True, True, True, True, True, True, + True, True, True, True, True, False, False, False, + False, False, False, False]]) + + assert np.allclose(evp.fields['reflectivity']['data'][0,10:30], evp_refl) + assert np.allclose(evp.fields['reflectivity']['data'][0,10:30].mask, evp_refl_mask) + + assert len(evp.range['data']) == 40 + assert np.all(evp.range['data'] == np.arange(125, 10000, 250)) + assert evp.azimuth['data'][0] == 0 + def test_quasi_vertical_profile(): test_radar = pyart.testing.make_target_radar() height = np.arange(0, 1000, 200) diff --git a/tests/util/test_circular_stats.py b/tests/util/test_circular_stats.py index 79063ec3848..6a8251f5f1a 100644 --- a/tests/util/test_circular_stats.py +++ b/tests/util/test_circular_stats.py @@ -6,6 +6,34 @@ from pyart.util import circular_stats +def test_compute_directional_stats(): + field_1d = np.ma.array([1,1,3,4,5]) + field_1d.mask = [True, False, False, False, False] + field = np.tile(field_1d, (10, 1)) + + mean, nvalid = circular_stats.compute_directional_stats(field, axis = 1) + median, nvalid = circular_stats.compute_directional_stats(field, axis = 1, + avg_type = 'median') + + assert np.all(mean == 3.25) + assert np.all(median == 3.5) + assert np.all(nvalid == 4) + + # Test with larger nb of nvalid_min + mean, nvalid = circular_stats.compute_directional_stats(field, axis = 1, + nvalid_min = 5) + assert np.all(mean.mask) + + # Test other direction + mean, nvalid = circular_stats.compute_directional_stats(field, axis = 0) + median, nvalid = circular_stats.compute_directional_stats(field, axis = 0, + avg_type = 'median') + + ref = np.ma.array([np.nan, 1, 3, 4, 5], mask = [1, 0 , 0, 0, 0]) + assert np.all(mean == ref) + assert np.all(median == ref) + assert np.all(nvalid == [0, 10, 10, 10, 10]) + def test_mean_of_two_angles(): mean = circular_stats.mean_of_two_angles(np.pi / 4 + 0.2, 7 * np.pi / 4) assert_almost_equal(mean, 0.10, 2) diff --git a/tests/util/test_radar_utils.py b/tests/util/test_radar_utils.py index f6dd40fdb1d..49d005eaf7a 100644 --- a/tests/util/test_radar_utils.py +++ b/tests/util/test_radar_utils.py @@ -104,6 +104,14 @@ def test_subset_radar(): assert radarcut.elevation["data"].max() <= ele_max assert list(radarcut.fields) == ["f1"] +def test_ma_broadcast_to(): + buf = np.ma.zeros(5) + buf.mask = [1,1,0,0,0] + buf_broad = pyart.util.radar_utils.ma_broadcast_to(buf, (10, 5)) + assert buf_broad.shape == (10,5) + assert buf_broad.mask.shape == (10,5) + expected_mask = np.tile(np.array([True, True, False, False, False]),[10,1]) + assert np.all(expected_mask == buf_broad.mask) # read in example file radar = pyart.io.read_nexrad_archive(pyart.testing.NEXRAD_ARCHIVE_MSG31_FILE) @@ -129,3 +137,4 @@ def test_image_mute_radar(): n_nonmutez = np.sum(~mute_radar.fields["nonmuted_reflectivity"]["data"].mask) assert n_mutez + n_nonmutez == n_rhohv + From d5b307231d0d41d3b7558147c033d5e4d72cb784 Mon Sep 17 00:00:00 2001 From: Daniel Wolfensberger Date: Tue, 27 Aug 2024 18:50:20 +0200 Subject: [PATCH 02/15] FIX: fixing est_rain_rate_hydro to match ARM original --- pyart/retrieve/qpe.py | 140 ++++++++++++++++++++++-------------------- 1 file changed, 72 insertions(+), 68 deletions(-) diff --git a/pyart/retrieve/qpe.py b/pyart/retrieve/qpe.py index 27bfd9ca5f2..48338791b9f 100644 --- a/pyart/retrieve/qpe.py +++ b/pyart/retrieve/qpe.py @@ -416,11 +416,23 @@ def est_rain_rate_za( return rain_main -def est_rain_rate_hydro(radar, alphazr=0.0376, betazr=0.6112, alphazs=0.1, - betazs=0.5, alphaa=None, betaa=None, mp_factor=0.6, - refl_field=None, a_field=None, hydro_field=None, - rr_field=None, master_field=None, thresh=None, - thresh_max=False): +def est_rain_rate_hydro( + radar, + alphazr=0.0376, + betazr=0.6112, + alphazs=0.1, + betazs=0.5, + alphaa=None, + betaa=None, + mp_factor=0.6, + refl_field=None, + a_field=None, + hydro_field=None, + rr_field=None, + main_field=None, + thresh=None, + thresh_max=False, +): """ Estimates rainfall rate using different relations between R and the polarimetric variables depending on the hydrometeor type. @@ -447,15 +459,15 @@ def est_rain_rate_hydro(radar, alphazr=0.0376, betazr=0.6112, alphazs=0.1, Name of the hydrometeor classification field to use. rr_field : str, optional Name of the rainfall rate field. - master_field : str, optional - Name of the field that is going to act as master. Has to be + main_field : str, optional + Name of the field that is going to act as main. Has to be either refl_field or kdp_field. Default is refl_field. thresh : float, optional - Value of the threshold that determines when to use the slave + Value of the threshold that determines when to use the secondary field. thresh_max : Bool, optional - If true the master field is used up to the thresh value maximum. - Otherwise the master field is not used below thresh value. + If true the main field is used up to the thresh value maximum. + Otherwise the main field is not used below thresh value. Returns ------- @@ -465,98 +477,90 @@ def est_rain_rate_hydro(radar, alphazr=0.0376, betazr=0.6112, alphazs=0.1, """ # parse the field parameters if refl_field is None: - refl_field = get_field_name('reflectivity') + refl_field = get_field_name("reflectivity") if a_field is None: - a_field = get_field_name('specific_attenuation') + a_field = get_field_name("specific_attenuation") if hydro_field is None: - hydro_field = get_field_name('radar_echo_classification') + hydro_field = get_field_name("radar_echo_classification") if rr_field is None: - rr_field = get_field_name('radar_estimated_rain_rate') - if master_field is None: - master_field = refl_field + rr_field = get_field_name("radar_estimated_rain_rate") # extract fields and parameters from radar if hydro_field in radar.fields: - hydroclass = radar.fields[hydro_field]['data'] + hydroclass = radar.fields[hydro_field]["data"] else: - raise KeyError('Field not available: ' + hydro_field) + raise KeyError("Field not available: " + hydro_field) # get the location of each hydrometeor class - is_ds = hydroclass == 2 - is_cr = hydroclass == 3 - is_lr = hydroclass == 4 - is_gr = hydroclass == 5 - is_rn = hydroclass == 6 - is_vi = hydroclass == 7 - is_ws = hydroclass == 8 - is_mh = hydroclass == 9 - is_ih = hydroclass == 10 + is_ds = hydroclass == 1 + is_cr = hydroclass == 2 + is_lr = hydroclass == 3 + is_gr = hydroclass == 4 + is_rn = hydroclass == 5 + is_vi = hydroclass == 6 + is_ws = hydroclass == 7 + is_mh = hydroclass == 8 + is_ih = hydroclass == 9 # compute z-r (in rain) z-r in snow and z-a relations rain_z = est_rain_rate_z( - radar, alpha=alphazr, beta=betazr, refl_field=refl_field, - rr_field=rr_field) + radar, alpha=alphazr, beta=betazr, refl_field=refl_field, rr_field=rr_field + ) snow_z = est_rain_rate_z( - radar, alpha=alphazs, beta=betazs, refl_field=refl_field, - rr_field=rr_field) + radar, alpha=alphazs, beta=betazs, refl_field=refl_field, rr_field=rr_field + ) rain_a = est_rain_rate_a( - radar, alpha=alphaa, beta=betaa, a_field=a_field, rr_field=rr_field) + radar, alpha=alphaa, beta=betaa, a_field=a_field, rr_field=rr_field + ) # initialize rainfall rate field - rr_data = np.ma.zeros(hydroclass.shape, dtype='float32') + rr_data = np.ma.zeros(hydroclass.shape, dtype="float32") rr_data[:] = np.ma.masked rr_data.set_fill_value(get_fillvalue()) # apply the relations for each type # solid phase - rr_data[is_ds] = snow_z['data'][is_ds] - rr_data[is_cr] = snow_z['data'][is_cr] - rr_data[is_vi] = snow_z['data'][is_vi] - rr_data[is_gr] = snow_z['data'][is_gr] - rr_data[is_ih] = snow_z['data'][is_ih] + rr_data[is_ds] = snow_z["data"][is_ds] + rr_data[is_cr] = snow_z["data"][is_cr] + rr_data[is_vi] = snow_z["data"][is_vi] + rr_data[is_gr] = snow_z["data"][is_gr] + rr_data[is_ih] = snow_z["data"][is_ih] # rain - if master_field == refl_field: - rain_master = rain_z - rain_slave = rain_a - if thresh is None: - thresh = 5. - thresh_max = True - elif master_field == a_field: - rain_master = rain_a - rain_slave = rain_z - if thresh is None: - thresh = 5. - thresh_max = False - elif master_field is None: - rain_master = rain_a - rain_slave = rain_z - thresh = 5. - thresh_max = False + if main_field == refl_field: + rain_main = rain_z + rain_secondary = rain_a + elif main_field == a_field: + rain_main = rain_a + rain_secondary = rain_z + elif main_field is None: + main_field = a_field + rain_main = rain_a + rain_secondary = rain_z else: - rain_master = rain_a - rain_slave = rain_z - thresh = 5. + main_field = a_field + rain_main = rain_a + rain_secondary = rain_z + thresh = 0.04 thresh_max = False - warn('Unknown master field. Using ' + a_field + ' with threshold ' + - str(thresh)) + warn("Unknown main field. Using " + a_field + " with threshold " + str(thresh)) if thresh_max: - is_slave = rain_master['data'] > thresh + is_secondary = rain_main["data"] > thresh else: - is_slave = rain_master['data'] < thresh + is_secondary = rain_main["data"] < thresh - rain_master['data'][is_slave] = rain_slave['data'][is_slave] + rain_main["data"][is_secondary] = rain_secondary["data"][is_secondary] - rr_data[is_lr] = rain_master['data'][is_lr] - rr_data[is_rn] = rain_master['data'][is_rn] + rr_data[is_lr] = rain_main["data"][is_lr] + rr_data[is_rn] = rain_main["data"][is_rn] # mixed phase - rr_data[is_ws] = mp_factor * rain_z['data'][is_ws] - rr_data[is_mh] = mp_factor * rain_z['data'][is_mh] + rr_data[is_ws] = mp_factor * rain_z["data"][is_ws] + rr_data[is_mh] = mp_factor * rain_z["data"][is_mh] rain = get_metadata(rr_field) - rain['data'] = rr_data + rain["data"] = rr_data return rain From a4f4fb075d65acffda265225fbe0ed62f20967fb Mon Sep 17 00:00:00 2001 From: Daniel Wolfensberger Date: Tue, 27 Aug 2024 19:01:32 +0200 Subject: [PATCH 03/15] ENH: black formatting --- examples/correct/plot_attenuation.py | 1 + examples/correct/plot_dealias.py | 1 + examples/io/plot_nexrad_data_aws.py | 1 + examples/io/plot_older_nexrad_data_aws.py | 1 + .../plot_compare_two_radars_gatemapper.py | 1 + .../mapping/plot_grid_single_sweep_ppi.py | 1 + .../mapping/plot_map_one_radar_to_grid.py | 1 + .../mapping/plot_map_two_radars_to_grid.py | 1 + examples/plotting/plot_choose_a_colormap.py | 1 + examples/plotting/plot_cross_section.py | 1 + examples/plotting/plot_modify_colorbar.py | 1 + .../plotting/plot_nexrad_multiple_moments.py | 1 + examples/plotting/plot_nexrad_reflectivity.py | 1 + examples/plotting/plot_ppi_cfradial.py | 1 + examples/plotting/plot_ppi_mdv.py | 1 + examples/plotting/plot_rhi_cfradial.py | 1 + .../plotting/plot_rhi_cfradial_singlescan.py | 1 + ..._rhi_contours_differential_reflectivity.py | 1 + examples/plotting/plot_rhi_data_overlay.py | 1 + examples/plotting/plot_rhi_mdv.py | 1 + examples/plotting/plot_rhi_two_panel.py | 1 + examples/plotting/plot_xsect.py | 1 + pyart/__check_build/__init__.py | 1 + pyart/core/transforms.py | 8 +- pyart/core/wind_profile.py | 1 - pyart/correct/phase_proc.py | 6 +- pyart/graph/_cm_colorblind.py | 1 + pyart/graph/gridmapdisplay_basemap.py | 1 - pyart/io/mdv_grid.py | 1 - pyart/io/nexrad_common.py | 1 + pyart/io/nexrad_level2.py | 1 - pyart/io/nexrad_level3.py | 2 +- pyart/map/gates_to_grid.py | 1 + pyart/retrieve/__init__.py | 4 +- pyart/retrieve/_echo_class_wt.py | 1 - pyart/retrieve/echo_class.py | 238 +++-- pyart/retrieve/qpe.py | 1 + pyart/retrieve/qvp.py | 953 +++++++++++------- pyart/util/__init__.py | 4 +- pyart/util/circular_stats.py | 6 +- pyart/util/radar_utils.py | 25 +- pyart/xradar/accessor.py | 1 - tests/bridge/test_wradlib_bridge.py | 1 - tests/core/test_transforms.py | 11 +- tests/graph/test_cm.py | 1 - tests/graph/test_cm_colorblind.py | 1 - tests/graph/test_gridmapdisplay.py | 1 + tests/graph/test_gridmapdisplay_basemap.py | 1 + tests/graph/test_radarmapdisplay.py | 1 + tests/graph/test_radarmapdisplay_basemap.py | 1 + tests/retrieve/test_echo_class.py | 9 +- tests/retrieve/test_qvp.py | 207 +++- tests/util/test_circular_stats.py | 22 +- tests/util/test_radar_utils.py | 11 +- 54 files changed, 947 insertions(+), 599 deletions(-) diff --git a/examples/correct/plot_attenuation.py b/examples/correct/plot_attenuation.py index 469b2d049fa..6bdbb3ee745 100644 --- a/examples/correct/plot_attenuation.py +++ b/examples/correct/plot_attenuation.py @@ -7,6 +7,7 @@ for a polarimetric radar using a Z-PHI method implemented in Py-ART. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/correct/plot_dealias.py b/examples/correct/plot_dealias.py index 8687b22c81b..e33a3acb1f8 100644 --- a/examples/correct/plot_dealias.py +++ b/examples/correct/plot_dealias.py @@ -6,6 +6,7 @@ In this example doppler velocities are dealiased using the ial condition of the dealiasing, using the region-based dealiasing algorithm in Py-ART. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/io/plot_nexrad_data_aws.py b/examples/io/plot_nexrad_data_aws.py index 7a0a055191c..a9eabc55412 100644 --- a/examples/io/plot_nexrad_data_aws.py +++ b/examples/io/plot_nexrad_data_aws.py @@ -7,6 +7,7 @@ and plot quick looks of the datasets. """ + print(__doc__) # Author: Max Grover (mgrover@anl.gov) diff --git a/examples/io/plot_older_nexrad_data_aws.py b/examples/io/plot_older_nexrad_data_aws.py index b50f8d81c64..e650da6fb3f 100644 --- a/examples/io/plot_older_nexrad_data_aws.py +++ b/examples/io/plot_older_nexrad_data_aws.py @@ -7,6 +7,7 @@ to 2008 that are missing some coordinate metadata. """ + print(__doc__) diff --git a/examples/mapping/plot_compare_two_radars_gatemapper.py b/examples/mapping/plot_compare_two_radars_gatemapper.py index fe57ea4557f..8b8274218ee 100644 --- a/examples/mapping/plot_compare_two_radars_gatemapper.py +++ b/examples/mapping/plot_compare_two_radars_gatemapper.py @@ -7,6 +7,7 @@ another radar in Antenna coordinates and compare the fields. """ + print(__doc__) # Author: Max Grover (mgrover@anl.gov) and Bobby Jackson (rjackson@anl.gov) diff --git a/examples/mapping/plot_grid_single_sweep_ppi.py b/examples/mapping/plot_grid_single_sweep_ppi.py index f7c5eca8ac6..beac728837d 100644 --- a/examples/mapping/plot_grid_single_sweep_ppi.py +++ b/examples/mapping/plot_grid_single_sweep_ppi.py @@ -13,6 +13,7 @@ negligible at low elevation angles common to PPI sweeps. """ + print(__doc__) # ===================== diff --git a/examples/mapping/plot_map_one_radar_to_grid.py b/examples/mapping/plot_map_one_radar_to_grid.py index fc2b9dfd561..9b5280ee7dd 100644 --- a/examples/mapping/plot_map_one_radar_to_grid.py +++ b/examples/mapping/plot_map_one_radar_to_grid.py @@ -7,6 +7,7 @@ Cartesian grid. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/mapping/plot_map_two_radars_to_grid.py b/examples/mapping/plot_map_two_radars_to_grid.py index 7298ea43568..b4b0c1cbfbb 100644 --- a/examples/mapping/plot_map_two_radars_to_grid.py +++ b/examples/mapping/plot_map_two_radars_to_grid.py @@ -7,6 +7,7 @@ coordinates to a Cartesian grid. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/plotting/plot_choose_a_colormap.py b/examples/plotting/plot_choose_a_colormap.py index acc1e9880df..c322a93442f 100644 --- a/examples/plotting/plot_choose_a_colormap.py +++ b/examples/plotting/plot_choose_a_colormap.py @@ -7,6 +7,7 @@ and how to add them to your own plots. """ + print(__doc__) # Author: Max Grover (mgrover@anl.gov) diff --git a/examples/plotting/plot_cross_section.py b/examples/plotting/plot_cross_section.py index aafa9e3b790..420d4b40153 100644 --- a/examples/plotting/plot_cross_section.py +++ b/examples/plotting/plot_cross_section.py @@ -7,6 +7,7 @@ of your radar grid using the GridMapDisplay """ + print(__doc__) # Author: Max Grover (mgrover@anl.gov) diff --git a/examples/plotting/plot_modify_colorbar.py b/examples/plotting/plot_modify_colorbar.py index 9451f897ca7..44495da9503 100644 --- a/examples/plotting/plot_modify_colorbar.py +++ b/examples/plotting/plot_modify_colorbar.py @@ -7,6 +7,7 @@ within a Py-ART display object. """ + print(__doc__) # Author: Joe O'Brien (obrienj@anl.gov) diff --git a/examples/plotting/plot_nexrad_multiple_moments.py b/examples/plotting/plot_nexrad_multiple_moments.py index b6c03fa9a26..7e6cc6041d7 100644 --- a/examples/plotting/plot_nexrad_multiple_moments.py +++ b/examples/plotting/plot_nexrad_multiple_moments.py @@ -7,6 +7,7 @@ NEXRAD Archive file. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/plotting/plot_nexrad_reflectivity.py b/examples/plotting/plot_nexrad_reflectivity.py index a3e73f0edca..49f12bccf39 100644 --- a/examples/plotting/plot_nexrad_reflectivity.py +++ b/examples/plotting/plot_nexrad_reflectivity.py @@ -7,6 +7,7 @@ NEXRAD file. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/plotting/plot_ppi_cfradial.py b/examples/plotting/plot_ppi_cfradial.py index 982e6602261..33f18ea8542 100644 --- a/examples/plotting/plot_ppi_cfradial.py +++ b/examples/plotting/plot_ppi_cfradial.py @@ -6,6 +6,7 @@ An example which creates a PPI plot of a Cfradial file. """ + print(__doc__) # Author: Max Grover (mgrover@anl.gov) diff --git a/examples/plotting/plot_ppi_mdv.py b/examples/plotting/plot_ppi_mdv.py index 706495bf3a6..2b1036e41f6 100644 --- a/examples/plotting/plot_ppi_mdv.py +++ b/examples/plotting/plot_ppi_mdv.py @@ -6,6 +6,7 @@ An example which creates a PPI plot of a MDV file using a RadarDisplay object. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/plotting/plot_rhi_cfradial.py b/examples/plotting/plot_rhi_cfradial.py index 5878a287064..fa85dd5fe09 100755 --- a/examples/plotting/plot_rhi_cfradial.py +++ b/examples/plotting/plot_rhi_cfradial.py @@ -7,6 +7,7 @@ a RadarDisplay object. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/plotting/plot_rhi_cfradial_singlescan.py b/examples/plotting/plot_rhi_cfradial_singlescan.py index 088b74d2c37..9f32eeedcf5 100755 --- a/examples/plotting/plot_rhi_cfradial_singlescan.py +++ b/examples/plotting/plot_rhi_cfradial_singlescan.py @@ -7,6 +7,7 @@ a RadarDisplay object. """ + print(__doc__) import matplotlib.pyplot as plt diff --git a/examples/plotting/plot_rhi_contours_differential_reflectivity.py b/examples/plotting/plot_rhi_contours_differential_reflectivity.py index bde4afbb58a..f64a9d0b8ca 100644 --- a/examples/plotting/plot_rhi_contours_differential_reflectivity.py +++ b/examples/plotting/plot_rhi_contours_differential_reflectivity.py @@ -7,6 +7,7 @@ and adding differnential Reflectivity contours from the same MDV file. """ + print(__doc__) # Author: Cory Weber (cweber@anl.gov) diff --git a/examples/plotting/plot_rhi_data_overlay.py b/examples/plotting/plot_rhi_data_overlay.py index 2e6a901e29d..c65b95f9843 100644 --- a/examples/plotting/plot_rhi_data_overlay.py +++ b/examples/plotting/plot_rhi_data_overlay.py @@ -7,6 +7,7 @@ and adding Reflectivity contours from the same MDV file. """ + print(__doc__) # Author: Cory Weber (cweber@anl.gov) diff --git a/examples/plotting/plot_rhi_mdv.py b/examples/plotting/plot_rhi_mdv.py index c5cfcfea67a..e6e5ca9cf4b 100644 --- a/examples/plotting/plot_rhi_mdv.py +++ b/examples/plotting/plot_rhi_mdv.py @@ -6,6 +6,7 @@ An example which creates a RHI plot of a MDV file using a RadarDisplay object. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/examples/plotting/plot_rhi_two_panel.py b/examples/plotting/plot_rhi_two_panel.py index 9d6b84ae890..9efb26a60b8 100755 --- a/examples/plotting/plot_rhi_two_panel.py +++ b/examples/plotting/plot_rhi_two_panel.py @@ -7,6 +7,7 @@ included in the two panels are reflectivity and doppler velocity. """ + print(__doc__) # Author: Max Grover (mgrover@anl.gov) diff --git a/examples/plotting/plot_xsect.py b/examples/plotting/plot_xsect.py index a64fcca70c4..5bc92beed34 100644 --- a/examples/plotting/plot_xsect.py +++ b/examples/plotting/plot_xsect.py @@ -7,6 +7,7 @@ of PPI scans and plots both cross sections. """ + print(__doc__) # Author: Jonathan J. Helmus (jhelmus@anl.gov) diff --git a/pyart/__check_build/__init__.py b/pyart/__check_build/__init__.py index 7c050370481..610b381934c 100644 --- a/pyart/__check_build/__init__.py +++ b/pyart/__check_build/__init__.py @@ -1,6 +1,7 @@ """ Module to give helpful messages to the user that did not compile Py-ART properly. """ + import os INPLACE_MSG = """ diff --git a/pyart/core/transforms.py b/pyart/core/transforms.py index a4c02f88429..d37c2b3528d 100644 --- a/pyart/core/transforms.py +++ b/pyart/core/transforms.py @@ -77,6 +77,7 @@ def antenna_to_cartesian(ranges, azimuths, elevations): y = s * np.cos(theta_a) return x, y, z + def cartesian_to_antenna(x, y, z): """ Returns antenna coordinates from Cartesian coordinates. @@ -98,13 +99,14 @@ def cartesian_to_antenna(x, y, z): Effective radius scale factor """ - ranges = np.sqrt(x ** 2. + y ** 2. + z ** 2.) - elevations = np.rad2deg(np.arctan(z / np.sqrt(x ** 2. + y ** 2.))) + ranges = np.sqrt(x**2.0 + y**2.0 + z**2.0) + elevations = np.rad2deg(np.arctan(z / np.sqrt(x**2.0 + y**2.0))) azimuths = np.rad2deg(np.arctan2(x, y)) # [-180, 180] - azimuths[azimuths < 0.] += 360. # [0, 360] + azimuths[azimuths < 0.0] += 360.0 # [0, 360] return ranges, azimuths, elevations + def antenna_vectors_to_cartesian(ranges, azimuths, elevations, edges=False): """ Calculate Cartesian coordinate for gates from antenna coordinate vectors. diff --git a/pyart/core/wind_profile.py b/pyart/core/wind_profile.py index 2a8520a79a0..2ebfad767cb 100644 --- a/pyart/core/wind_profile.py +++ b/pyart/core/wind_profile.py @@ -3,7 +3,6 @@ """ - import numpy as np diff --git a/pyart/correct/phase_proc.py b/pyart/correct/phase_proc.py index 27b6901f1c9..2d76473035d 100644 --- a/pyart/correct/phase_proc.py +++ b/pyart/correct/phase_proc.py @@ -108,11 +108,7 @@ def fzl_index(fzl, ranges, elevation, radar_height): p_r = 4.0 * Re / 3.0 z = ( radar_height - + ( - ranges**2 - + p_r**2 - + 2.0 * ranges * p_r * np.sin(elevation * np.pi / 180.0) - ) + + (ranges**2 + p_r**2 + 2.0 * ranges * p_r * np.sin(elevation * np.pi / 180.0)) ** 0.5 - p_r ) diff --git a/pyart/graph/_cm_colorblind.py b/pyart/graph/_cm_colorblind.py index 8b078ea58e9..164fcd02dec 100644 --- a/pyart/graph/_cm_colorblind.py +++ b/pyart/graph/_cm_colorblind.py @@ -2,6 +2,7 @@ Data for colorblind friendly radar colormaps """ + import os import numpy as np diff --git a/pyart/graph/gridmapdisplay_basemap.py b/pyart/graph/gridmapdisplay_basemap.py index 1b6f914aaee..997de186ea3 100644 --- a/pyart/graph/gridmapdisplay_basemap.py +++ b/pyart/graph/gridmapdisplay_basemap.py @@ -3,7 +3,6 @@ """ - import matplotlib.pyplot as plt import numpy as np diff --git a/pyart/io/mdv_grid.py b/pyart/io/mdv_grid.py index ebcedf3bb87..e8619665001 100644 --- a/pyart/io/mdv_grid.py +++ b/pyart/io/mdv_grid.py @@ -3,7 +3,6 @@ """ - import datetime import warnings diff --git a/pyart/io/nexrad_common.py b/pyart/io/nexrad_common.py index e34c27751e2..13f3c6acd8a 100644 --- a/pyart/io/nexrad_common.py +++ b/pyart/io/nexrad_common.py @@ -2,6 +2,7 @@ Data and functions common to all types of NEXRAD files. """ + # The functions in this module are intended to be used in other # nexrad related modules. The functions are not and should not be # exported into the pyart.io namespace. diff --git a/pyart/io/nexrad_level2.py b/pyart/io/nexrad_level2.py index 38a3e88530b..0facc446dbd 100644 --- a/pyart/io/nexrad_level2.py +++ b/pyart/io/nexrad_level2.py @@ -2,7 +2,6 @@ Functions for reading NEXRAD level 2 files. """ - # This file is part of the Py-ART, the Python ARM Radar Toolkit # https://github.com/ARM-DOE/pyart diff --git a/pyart/io/nexrad_level3.py b/pyart/io/nexrad_level3.py index 0b2ed137f47..72605fb3a6e 100644 --- a/pyart/io/nexrad_level3.py +++ b/pyart/io/nexrad_level3.py @@ -761,7 +761,7 @@ def _int16_to_float16(val): ("block_length", INT4), # Length of block in bytes ("layers", INT2), # Number of data layers ("layer_divider", INT2), # Delineate data layers, -1 - ("layer_length", INT4) # Length of data layer in bytes + ("layer_length", INT4), # Length of data layer in bytes # Display data packets ) diff --git a/pyart/map/gates_to_grid.py b/pyart/map/gates_to_grid.py index bcc470d55a9..95a6e888248 100644 --- a/pyart/map/gates_to_grid.py +++ b/pyart/map/gates_to_grid.py @@ -2,6 +2,7 @@ Generate a Cartesian grid by mapping from radar gates onto the grid. """ + import gc import warnings diff --git a/pyart/retrieve/__init__.py b/pyart/retrieve/__init__.py index 60893b01f29..563ced9a2cb 100644 --- a/pyart/retrieve/__init__.py +++ b/pyart/retrieve/__init__.py @@ -21,8 +21,8 @@ from .qpe import est_rain_rate_zkdp # noqa from .qpe import est_rain_rate_zpoly # noqa from .qpe import ZtoR # noqa -from .qvp import quasi_vertical_profile, compute_qvp, compute_rqvp #noqa -from .qvp import compute_evp, compute_svp, compute_vp, compute_ts_along_coord #noqa +from .qvp import quasi_vertical_profile, compute_qvp, compute_rqvp # noqa +from .qvp import compute_evp, compute_svp, compute_vp, compute_ts_along_coord # noqa from .simple_moment_calculations import calculate_snr_from_reflectivity # noqa from .simple_moment_calculations import calculate_velocity_texture # noqa from .simple_moment_calculations import compute_cdr # noqa diff --git a/pyart/retrieve/_echo_class_wt.py b/pyart/retrieve/_echo_class_wt.py index e9d38365e34..1e1c30615f7 100644 --- a/pyart/retrieve/_echo_class_wt.py +++ b/pyart/retrieve/_echo_class_wt.py @@ -13,7 +13,6 @@ atwt2d """ - import numpy as np diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py index 815e5149f08..5cb61e73024 100644 --- a/pyart/retrieve/echo_class.py +++ b/pyart/retrieve/echo_class.py @@ -601,19 +601,27 @@ def feature_detection( return feature_dict -def hydroclass_semisupervised(radar, - hydro_names=('AG', 'CR', 'LR', 'RP', 'RN', 'VI', - 'WS', 'MH', 'IH/HDG'), - var_names=('Zh', 'ZDR', 'KDP', 'RhoHV', - 'relH'), - mass_centers=None, - weights=np.array([1., 1., 1., 0.75, 0.5]), - value=50., lapse_rate=-6.5, refl_field=None, - zdr_field=None, rhv_field=None, kdp_field=None, - temp_field=None, iso0_field=None, - hydro_field=None, entropy_field=None, - temp_ref='temperature', compute_entropy=False, - output_distances=False, vectorize=False): +def hydroclass_semisupervised( + radar, + hydro_names=("AG", "CR", "LR", "RP", "RN", "VI", "WS", "MH", "IH/HDG"), + var_names=("Zh", "ZDR", "KDP", "RhoHV", "relH"), + mass_centers=None, + weights=np.array([1.0, 1.0, 1.0, 0.75, 0.5]), + value=50.0, + lapse_rate=-6.5, + refl_field=None, + zdr_field=None, + rhv_field=None, + kdp_field=None, + temp_field=None, + iso0_field=None, + hydro_field=None, + entropy_field=None, + temp_ref="temperature", + compute_entropy=False, + output_distances=False, + vectorize=False, +): """ Classifies precipitation echoes into hydrometeor types. @@ -665,7 +673,7 @@ def hydroclass_semisupervised(radar, ------- fields_dict : dict Dictionary containing the retrieved fields - + The output directionary field_dict has the following keys: hydro : dict @@ -680,14 +688,14 @@ def hydroclass_semisupervised(radar, - 7: Wet snow - 8: Melting hail - 9: Dry hail or high-density graupel - + if compute_entropy is True: entropy : dict Shannon entropy of the hydrometeor demixing propX : dict - Proportion of a given hydrometeor class in the polarimetric + Proportion of a given hydrometeor class in the polarimetric decomposition of a radar volume References @@ -697,7 +705,7 @@ def hydroclass_semisupervised(radar, of polarimetric radar measurements: a semi-supervised approach, Atmos. Meas. Tech., 9, 4425-4445, doi:10.5194/amt-9-4425-2016, 2016 - Besic, N., Gehring, J., Praz, C., Figueras i Ventura, J., Grazioli, J., Gabella, + Besic, N., Gehring, J., Praz, C., Figueras i Ventura, J., Grazioli, J., Gabella, M., Germann, U., and Berne, A.: Unraveling hydrometeor mixtures in polarimetric radar measurements, Atmos. Meas. Tech., 11, 4847–4866, https://doi.org/10.5194/amt-11-4847-2018, 2018. @@ -731,74 +739,76 @@ def hydroclass_semisupervised(radar, "data": [9.2e9] } """ - + # select the centroids as a function of frequency band if mass_centers is None: # assign coefficients according to radar frequency - if 'frequency' in radar.instrument_parameters: + if "frequency" in radar.instrument_parameters: mass_centers = _get_mass_centers( - radar.instrument_parameters['frequency']['data'][0]) + radar.instrument_parameters["frequency"]["data"][0] + ) else: - mass_centers = _mass_centers_table()['C'] - warn('Radar frequency unknown. ' + - 'Default coefficients for C band will be applied') + mass_centers = _mass_centers_table()["C"] + warn( + "Radar frequency unknown. " + + "Default coefficients for C band will be applied" + ) if hydro_field is None: - hydro_field = get_field_name('radar_echo_classification') + hydro_field = get_field_name("radar_echo_classification") if compute_entropy: if entropy_field is None: - entropy_field = get_field_name('hydroclass_entropy') + entropy_field = get_field_name("hydroclass_entropy") # Get the data fields fields_dict = {} for var_name in var_names: - if var_name == 'Zh': + if var_name == "Zh": if refl_field is None: - refl_field = get_field_name('reflectivity') + refl_field = get_field_name("reflectivity") radar.check_field_exists(refl_field) - fields_dict.update({ - var_name: deepcopy(radar.fields[refl_field]['data'])}) - elif var_name == 'ZDR': + fields_dict.update({var_name: deepcopy(radar.fields[refl_field]["data"])}) + elif var_name == "ZDR": if zdr_field is None: - zdr_field = get_field_name('differential_reflectivity') + zdr_field = get_field_name("differential_reflectivity") radar.check_field_exists(zdr_field) - fields_dict.update({ - var_name: deepcopy(radar.fields[zdr_field]['data'])}) - elif var_name == 'KDP': + fields_dict.update({var_name: deepcopy(radar.fields[zdr_field]["data"])}) + elif var_name == "KDP": if kdp_field is None: - kdp_field = get_field_name('specific_differential_phase') + kdp_field = get_field_name("specific_differential_phase") radar.check_field_exists(kdp_field) - fields_dict.update({ - var_name: deepcopy(radar.fields[kdp_field]['data'])}) - elif var_name == 'RhoHV': + fields_dict.update({var_name: deepcopy(radar.fields[kdp_field]["data"])}) + elif var_name == "RhoHV": if rhv_field is None: - rhv_field = get_field_name('cross_correlation_ratio') + rhv_field = get_field_name("cross_correlation_ratio") radar.check_field_exists(rhv_field) - fields_dict.update({ - var_name: deepcopy(radar.fields[rhv_field]['data'])}) - elif var_name == 'relH': - if temp_ref == 'temperature': + fields_dict.update({var_name: deepcopy(radar.fields[rhv_field]["data"])}) + elif var_name == "relH": + if temp_ref == "temperature": if temp_field is None: - temp_field = get_field_name('temperature') + temp_field = get_field_name("temperature") radar.check_field_exists(temp_field) # convert temp in relative height respect to iso0 - temp = deepcopy(radar.fields[temp_field]['data']) - fields_dict.update({var_name: temp * (1000. / lapse_rate)}) + temp = deepcopy(radar.fields[temp_field]["data"]) + fields_dict.update({var_name: temp * (1000.0 / lapse_rate)}) else: if iso0_field is None: - iso0_field = get_field_name('height_over_iso0') + iso0_field = get_field_name("height_over_iso0") radar.check_field_exists(iso0_field) - fields_dict.update({ - var_name: deepcopy(radar.fields[iso0_field]['data'])}) + fields_dict.update( + {var_name: deepcopy(radar.fields[iso0_field]["data"])} + ) else: raise ValueError( - 'Variable ' + var_name + ' unknown. ' - + 'Valid variable names for hydrometeor classification are: ' - + 'relH, Zh, ZDR, KDP and RhoHV') + "Variable " + + var_name + + " unknown. " + + "Valid variable names for hydrometeor classification are: " + + "relH, Zh, ZDR, KDP and RhoHV" + ) # standardize data and centroids - mc_std = np.empty( - np.shape(mass_centers), dtype=fields_dict[var_names[0]].dtype) + mc_std = np.empty(np.shape(mass_centers), dtype=fields_dict[var_names[0]].dtype) for i, var_name in enumerate(var_names): mc_std[:, i] = _standardize(mass_centers[:, i], var_name) fields_dict[var_name] = _standardize(fields_dict[var_name], var_name) @@ -806,46 +816,43 @@ def hydroclass_semisupervised(radar, # if entropy has to be computed get transformation parameters t_vals = None if compute_entropy: - t_vals = _compute_coeff_transform( - mc_std, weights=weights, value=value) + t_vals = _compute_coeff_transform(mc_std, weights=weights, value=value) # assign to class if vectorize: hydroclass_data, entropy_data, prop_data = _assign_to_class_scan( - fields_dict, mc_std, var_names=var_names, weights=weights, - t_vals=t_vals) + fields_dict, mc_std, var_names=var_names, weights=weights, t_vals=t_vals + ) else: hydroclass_data, entropy_data, prop_data = _assign_to_class( - fields_dict, mc_std, weights=weights, t_vals=t_vals) + fields_dict, mc_std, weights=weights, t_vals=t_vals + ) # prepare output fields fields_dict = dict() hydro = get_metadata(hydro_field) - hydro['data'] = hydroclass_data - hydro.update({'_FillValue': 0}) - labels = ['NC'] + hydro["data"] = hydroclass_data + hydro.update({"_FillValue": 0}) + labels = ["NC"] ticks = [1] boundaries = [0.5, 1.5] for i, hydro_name in enumerate(hydro_names): labels.append(hydro_name) ticks.append(i + 2) boundaries.append(i + 2.5) - hydro.update({ - 'labels': labels, - 'ticks': ticks, - 'boundaries': boundaries}) - fields_dict.update({'hydro': hydro}) + hydro.update({"labels": labels, "ticks": ticks, "boundaries": boundaries}) + fields_dict.update({"hydro": hydro}) if compute_entropy: entropy = get_metadata(entropy_field) - entropy['data'] = entropy_data - fields_dict.update({'entropy': entropy}) + entropy["data"] = entropy_data + fields_dict.update({"entropy": entropy}) if output_distances: for field_name in hydro_names: - field_name = 'proportion_' + field_name + field_name = "proportion_" + field_name prop = get_metadata(field_name) - prop['data'] = prop_data[:, :, i] + prop["data"] = prop_data[:, :, i] fields_dict.update({field_name: prop}) return fields_dict @@ -900,10 +907,14 @@ def _standardize(data, field_name, mx=None, mn=None): return field_std -def _assign_to_class(fields_dict, mass_centers, - var_names=('Zh', 'ZDR', 'KDP', 'RhoHV', 'relH'), - weights=np.array([1., 1., 1., 0.75, 0.5]), - t_vals=None): + +def _assign_to_class( + fields_dict, + mass_centers, + var_names=("Zh", "ZDR", "KDP", "RhoHV", "relH"), + weights=np.array([1.0, 1.0, 1.0, 0.75, 0.5]), + t_vals=None, +): """ Assigns an hydrometeor class to a radar range bin computing the distance between the radar variables an a centroid. @@ -953,7 +964,8 @@ def _assign_to_class(fields_dict, mass_centers, data.append(fields_dict[var_name][ray, :]) data = np.ma.array(data, dtype=dtype) weights_mat = np.broadcast_to( - weights.reshape(nvariables, 1), (nvariables, nbins)) + weights.reshape(nvariables, 1), (nvariables, nbins) + ) dist = np.ma.zeros((nclasses, nbins), dtype=dtype) # compute distance: masked entries will not contribute to the distance @@ -961,9 +973,11 @@ def _assign_to_class(fields_dict, mass_centers, for i in range(nclasses): centroids_class = mass_centers[i, :] centroids_class = np.broadcast_to( - centroids_class.reshape(nvariables, 1), (nvariables, nbins)) - dist_ray = np.ma.sqrt(np.ma.sum( - ((centroids_class - data)**2.) * weights_mat, axis=0)) + centroids_class.reshape(nvariables, 1), (nvariables, nbins) + ) + dist_ray = np.ma.sqrt( + np.ma.sum(((centroids_class - data) ** 2.0) * weights_mat, axis=0) + ) dist_ray[mask] = np.ma.masked dist[i, :] = dist_ray @@ -978,34 +992,36 @@ def _assign_to_class(fields_dict, mass_centers, # Transform the distance using the coefficient of the dominant class t_vals_ray = np.ma.masked_where(mask, t_vals[class_vec[0, :]]) - t_vals_ray = ma_broadcast_to( - t_vals_ray.reshape(1, nbins), (nclasses, nbins)) + t_vals_ray = ma_broadcast_to(t_vals_ray.reshape(1, nbins), (nclasses, nbins)) t_dist_ray = np.ma.exp(-t_vals_ray * dist) # set transformed distances to a value between 0 and 1 dist_total = np.ma.sum(t_dist_ray, axis=0) - dist_total = ma_broadcast_to( - dist_total.reshape(1, nbins), (nclasses, nbins)) + dist_total = ma_broadcast_to(dist_total.reshape(1, nbins), (nclasses, nbins)) t_dist_ray /= dist_total # Compute entropy entropy_ray = -np.ma.sum( - t_dist_ray * np.ma.log(t_dist_ray) / np.ma.log(nclasses), axis=0) + t_dist_ray * np.ma.log(t_dist_ray) / np.ma.log(nclasses), axis=0 + ) entropy_ray[mask] = np.ma.masked entropy[ray, :] = entropy_ray t_dist[ray, :, :] = np.ma.transpose(t_dist_ray) if t_vals is not None: - t_dist *= 100. + t_dist *= 100.0 return hydroclass, entropy, t_dist -def _assign_to_class_scan(fields_dict, mass_centers, - var_names=('Zh', 'ZDR', 'KDP', 'RhoHV', 'relH'), - weights=np.array([1., 1., 1., 0.75, 0.5]), - t_vals=None): +def _assign_to_class_scan( + fields_dict, + mass_centers, + var_names=("Zh", "ZDR", "KDP", "RhoHV", "relH"), + weights=np.array([1.0, 1.0, 1.0, 0.75, 0.5]), + t_vals=None, +): """ assigns an hydrometeor class to a radar range bin computing the distance between the radar variables an a centroid. @@ -1049,7 +1065,8 @@ def _assign_to_class_scan(fields_dict, mass_centers, data.append(fields_dict[var_name]) data = np.ma.array(data, dtype=dtype) weights_mat = np.broadcast_to( - weights.reshape(nvariables, 1, 1), (nvariables, nrays, nbins)) + weights.reshape(nvariables, 1, 1), (nvariables, nrays, nbins) + ) # compute distance: masked entries will not contribute to the distance mask = np.ma.getmaskarray(fields_dict[var_names[0]]) @@ -1059,10 +1076,11 @@ def _assign_to_class_scan(fields_dict, mass_centers, for i in range(nclasses): centroids_class = mass_centers[i, :] centroids_class = np.broadcast_to( - centroids_class.reshape(nvariables, 1, 1), - (nvariables, nrays, nbins)) - dist_aux = np.ma.sqrt(np.ma.sum( - ((centroids_class - data)**2.) * weights_mat, axis=0)) + centroids_class.reshape(nvariables, 1, 1), (nvariables, nrays, nbins) + ) + dist_aux = np.ma.sqrt( + np.ma.sum(((centroids_class - data) ** 2.0) * weights_mat, axis=0) + ) dist_aux[mask] = np.ma.masked dist[:, :, i] = dist_aux @@ -1078,26 +1096,28 @@ def _assign_to_class_scan(fields_dict, mass_centers, # Transform the distance using the coefficient of the dominant class t_vals_aux = np.ma.masked_where(mask, t_vals[class_vec[:, :, 0]]) t_vals_aux = ma_broadcast_to( - t_vals_aux.reshape(nrays, nbins, 1), (nrays, nbins, nclasses)) + t_vals_aux.reshape(nrays, nbins, 1), (nrays, nbins, nclasses) + ) t_dist = np.ma.exp(-t_vals_aux * dist) del t_vals_aux # set distance to a value between 0 and 1 dist_total = np.ma.sum(t_dist, axis=-1) dist_total = ma_broadcast_to( - dist_total.reshape(nrays, nbins, 1), (nrays, nbins, nclasses)) + dist_total.reshape(nrays, nbins, 1), (nrays, nbins, nclasses) + ) t_dist /= dist_total del dist_total # compute entroy - entropy = -np.ma.sum( - t_dist * np.ma.log(t_dist) / np.ma.log(nclasses), axis=-1) + entropy = -np.ma.sum(t_dist * np.ma.log(t_dist) / np.ma.log(nclasses), axis=-1) entropy[mask] = np.ma.masked - t_dist *= 100. + t_dist *= 100.0 return hydroclass, entropy, t_dist + def _get_mass_centers(freq): """ Get mass centers for a particular frequency. @@ -1400,9 +1420,10 @@ def conv_strat_raut( return reclass_dict -def _compute_coeff_transform(mass_centers, - weights=np.array([1., 1., 1., 0.75, 0.5]), - value=50.): + +def _compute_coeff_transform( + mass_centers, weights=np.array([1.0, 1.0, 1.0, 0.75, 0.5]), value=50.0 +): """ get the transformation coefficients @@ -1426,13 +1447,18 @@ def _compute_coeff_transform(mass_centers, t_vals = np.empty((nclasses, nclasses), dtype=mass_centers.dtype) for i in range(nclasses): weights_mat = np.broadcast_to( - weights.reshape(1, nvariables), (nclasses, nvariables)) + weights.reshape(1, nvariables), (nclasses, nvariables) + ) centroids_class = mass_centers[i, :] centroids_class = np.broadcast_to( - centroids_class.reshape(1, nvariables), (nclasses, nvariables)) + centroids_class.reshape(1, nvariables), (nclasses, nvariables) + ) t_vals[i, :] = np.sqrt( - np.sum(weights_mat * np.power( - np.abs(centroids_class - mass_centers), 2.), axis=1)) + np.sum( + weights_mat * np.power(np.abs(centroids_class - mass_centers), 2.0), + axis=1, + ) + ) # pick the second lowest value (the first is 0) t_vals = np.sort(t_vals, axis=-1)[:, 1] diff --git a/pyart/retrieve/qpe.py b/pyart/retrieve/qpe.py index 48338791b9f..c657228d826 100644 --- a/pyart/retrieve/qpe.py +++ b/pyart/retrieve/qpe.py @@ -416,6 +416,7 @@ def est_rain_rate_za( return rain_main + def est_rain_rate_hydro( radar, alphazr=0.0376, diff --git a/pyart/retrieve/qvp.py b/pyart/retrieve/qvp.py index 562f1236eb6..18ef9e52e12 100644 --- a/pyart/retrieve/qvp.py +++ b/pyart/retrieve/qvp.py @@ -41,8 +41,7 @@ from ..util.xsect import cross_section_ppi, cross_section_rhi -def quasi_vertical_profile(radar, desired_angle=None, fields=None, - gatefilter=None): +def quasi_vertical_profile(radar, desired_angle=None, fields=None, gatefilter=None): """ Quasi Vertical Profile. @@ -94,12 +93,12 @@ def quasi_vertical_profile(radar, desired_angle=None, fields=None, # radar angle if desired_angle is None: desired_angle = 20.0 - index = abs(radar.fixed_angle['data'] - desired_angle).argmin() + 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]) + 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 @@ -113,7 +112,8 @@ def quasi_vertical_profile(radar, desired_angle=None, fields=None, 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) + 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) @@ -126,7 +126,8 @@ def quasi_vertical_profile(radar, desired_angle=None, fields=None, 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) + 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) @@ -134,16 +135,27 @@ def quasi_vertical_profile(radar, desired_angle=None, fields=None, 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}) + 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 -def compute_qvp(radar, field_names, ref_time=None, angle=0., ang_tol=1., - hmax=10000., hres=50., avg_type='mean', nvalid_min=30, - interp_kind='none', qvp=None): +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. @@ -196,61 +208,86 @@ def compute_qvp(radar, field_names, ref_time=None, angle=0., ang_tol=1., Radar Data. JTECH vol. 33 pp 551-562 """ - if avg_type not in ('mean', 'median'): - warn('Unsuported statistics ' + avg_type) + 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': + if radar_aux.scan_type == "rhi": radar_aux = cross_section_rhi(radar_aux, [angle], el_tol=ang_tol) - elif radar_aux.scan_type == 'ppi': + elif radar_aux.scan_type == "ppi": radar_aux = radar_aux.extract_sweeps([int(angle)]) else: - warn('Error: unsupported scan type.') + 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) + 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]) + 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') + 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) + 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) + 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) + 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))) + qvp.fields[field_name]["data"] = np.ma.concatenate( + (qvp.fields[field_name]["data"], qvp_data.reshape(1, qvp.ngates)) + ) return qvp -def compute_rqvp(radar, field_names, ref_time=None, hmax=10000., hres=2., - avg_type='mean', nvalid_min=30, interp_kind='nearest', - rmax=50000., weight_power=2., qvp=None): +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. @@ -302,81 +339,96 @@ def compute_rqvp(radar, field_names, ref_time=None, hmax=10000., hres=2., Transitions. Weather and Forecasting vol. 32 pp 2065-2082 """ - if avg_type not in ('mean', 'median'): - warn('Unsuported statistics ' + avg_type) + 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': + 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 = 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.') + 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) + 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.) + qvp, ref_time, qvp.longitude["data"][0], qvp.latitude["data"][0], elev=90.0 + ) - rmax_km = rmax / 1000. + 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))}) + 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, :] + height = radar_aux.gate_altitude["data"][0, :] # compute ground range [Km] - grng = np.sqrt( - np.power(radar_aux.gate_x['data'][0, :], 2.) + - np.power(radar_aux.gate_y['data'][0, :], 2.)) / 1000. + 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']) + 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') + 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) + 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) + values, height, qvp.range["data"], interp_kind=interp_kind + ) # Compute weight - weight = np.ma.abs(grng_interp - (rmax_km - 1.)) - weight[grng_interp <= rmax_km - 1.] = 1. / np.power( - weight[grng_interp <= rmax_km - 1.], 0.) + 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. + weight[grng_interp > rmax_km - 1.0] = 0.0 else: - weight[grng_interp > rmax_km - 1.] = 1. / np.power( - weight[grng_interp > rmax_km - 1.], weight_power) + 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: @@ -385,25 +437,37 @@ def compute_rqvp(radar, field_names, ref_time=None, hmax=10000., hres=2., 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)) + 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) + 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))) + qvp.fields[field_name]["data"] = np.ma.concatenate( + (qvp.fields[field_name]["data"], qvp_data.reshape(1, qvp.ngates)) + ) return qvp -def compute_evp(radar, field_names, lon, lat, ref_time=None, - latlon_tol=0.0005, delta_rng=15000., delta_azi=10, - hmax=10000., hres=250., avg_type='mean', nvalid_min=1, - interp_kind='none', qvp=None): +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. @@ -457,33 +521,37 @@ def compute_evp(radar, field_names, lon, lat, ref_time=None, Meteorologische Zeitschrift vol. 26 pp 207-222 """ - if avg_type not in ('mean', 'median'): - warn('Unsuported statistics ' + avg_type) + 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': + 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 = 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.') + 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) + 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.) + qvp = _update_qvp_metadata(qvp, ref_time, lon, lat, elev=90.0) height = dict() values = dict() @@ -497,53 +565,73 @@ def compute_evp(radar, field_names, lon, lat, ref_time=None, # find nearest gate to lat lon point ind_ray, _, azi, rng = find_nearest_gate( - radar_aux, lat, lon, latlon_tol=latlon_tol) + 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) + 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') + 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]) + 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 = 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) + 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) + 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) + 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))) + qvp.fields[field_name]["data"] = np.ma.concatenate( + (qvp.fields[field_name]["data"], qvp_data.reshape(1, qvp.ngates)) + ) return qvp -def compute_svp(radar, field_names, lon, lat, angle, ref_time=None, - ang_tol=1., latlon_tol=0.0005, delta_rng=15000., delta_azi=10, - hmax=10000., hres=250., avg_type='mean', nvalid_min=1, - interp_kind='none', qvp=None): +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. @@ -603,87 +691,100 @@ def compute_svp(radar, field_names, lon, lat, angle, ref_time=None, in Central Oklahoma. JTECH vol. 56 pp 1345-1363 """ - if avg_type not in ('mean', 'median'): - warn('Unsuported statistics ' + avg_type) + 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': + 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.') + 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) + 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]) + 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) + 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') + 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) + 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))) + 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] + 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') + 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 = 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) + 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) + 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) + 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))) + qvp.fields[field_name]["data"] = np.ma.concatenate( + (qvp.fields[field_name]["data"], qvp_data.reshape(1, qvp.ngates)) + ) return qvp -def compute_vp(radar, field_names, lon, lat, ref_time=None, - latlon_tol=0.0005, hmax=10000., hres=50., - interp_kind='none', qvp=None): +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. @@ -726,25 +827,29 @@ def compute_vp(radar, field_names, lon, lat, ref_time=None, """ radar_aux = deepcopy(radar) # transform radar into ppi over the required elevation - if radar_aux.scan_type == 'rhi': + 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 = 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.') + 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) + 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.) + qvp = _update_qvp_metadata(qvp, ref_time, lon, lat, elev=90.0) height = dict() values = dict() @@ -757,45 +862,59 @@ def compute_vp(radar, field_names, lon, lat, ref_time=None, # find nearest gate to lat lon point ind_ray, ind_rng, _, _ = find_nearest_gate( - radar_aux, lat, lon, latlon_tol=latlon_tol) + 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') + 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]) + 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]) + 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) + 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) + 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))) + qvp.fields[field_name]["data"] = np.ma.concatenate( + (qvp.fields[field_name]["data"], qvp_data.reshape(1, qvp.ngates)) + ) return qvp -def compute_ts_along_coord(radar, field_names, mode='ALONG_AZI', - fixed_range=None, fixed_azimuth=None, - fixed_elevation=None, ang_tol=1., rng_tol=50., - value_start=None, value_stop=None, ref_time=None, - acoord=None): +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 @@ -831,38 +950,38 @@ def compute_ts_along_coord(radar, field_names, mode='ALONG_AZI', The computed data along a coordinate """ - if mode == 'ALONG_RNG': + if mode == "ALONG_RNG": if value_start is None: - value_start = 0. + value_start = 0.0 if value_stop is None: - value_stop = radar.range['data'][-1] + value_stop = radar.range["data"][-1] if fixed_azimuth is None: - fixed_azimuth = 0. + fixed_azimuth = 0.0 if fixed_elevation is None: - fixed_elevation = 0. - elif mode == 'ALONG_AZI': + fixed_elevation = 0.0 + elif mode == "ALONG_AZI": if value_start is None: - value_start = np.min(radar.azimuth['data']) + value_start = np.min(radar.azimuth["data"]) if value_stop is None: - value_stop = np.max(radar.azimuth['data']) + value_stop = np.max(radar.azimuth["data"]) if fixed_range is None: - fixed_range = 0. + fixed_range = 0.0 if fixed_elevation is None: - fixed_elevation = 0. - elif mode == 'ALONG_ELE': + fixed_elevation = 0.0 + elif mode == "ALONG_ELE": if value_start is None: - value_start = np.min(radar.elevation['data']) + value_start = np.min(radar.elevation["data"]) if value_stop is None: - value_stop = np.max(radar.elevation['data']) + value_stop = np.max(radar.elevation["data"]) if fixed_range is None: - fixed_range = 0. + fixed_range = 0.0 if fixed_azimuth is None: - fixed_azimuth = 0. + fixed_azimuth = 0.0 else: - warn('Unknown time series of coordinate mode ' + mode) + warn("Unknown time series of coordinate mode " + mode) return None - if mode == 'ALONG_RNG': + if mode == "ALONG_RNG": # rng_values : range # fixed_angle: elevation # elevation: elevation @@ -870,18 +989,24 @@ def compute_ts_along_coord(radar, field_names, mode='ALONG_AZI', 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) + 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') + 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': + elif mode == "ALONG_AZI": # rng_values : azimuth # fixed_angle : elevation # elevation : elevation @@ -889,12 +1014,18 @@ def compute_ts_along_coord(radar, field_names, mode='ALONG_AZI', 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) + 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') + warn("No data found") return None vals_dict.update({field_name: vals}) fixed_angle = fixed_elevation @@ -909,12 +1040,18 @@ def compute_ts_along_coord(radar, field_names, mode='ALONG_AZI', 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) + 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') + warn("No data found") return None vals_dict.update({field_name: vals}) fixed_angle = fixed_azimuth @@ -924,34 +1061,38 @@ def compute_ts_along_coord(radar, field_names, mode='ALONG_AZI', if acoord is None: acoord = _create_along_coord_object( - radar, field_names, rng_values, fixed_angle, mode, - start_time=ref_time) + 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') + 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) + 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) + 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))) + 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.): +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 @@ -978,33 +1119,41 @@ def project_to_vertical(data_in, data_height, grid_height, interp_kind='none', data_out = np.ma.masked_all(grid_height.size) return data_out - if interp_kind == 'none': + 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.) + 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': + 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_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) + 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.): +def find_rng_index(rng_vec, rng, rng_tol=0.0): """ Find the range index corresponding to a particular range @@ -1047,10 +1196,9 @@ def get_target_elevations(radar_in): 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]) + 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 @@ -1078,29 +1226,43 @@ def find_nearest_gate(radar, lat, lon, latlon_tol=0.0005): """ # find gates close to lat lon point - inds_ray_aux, inds_rng_aux = np.where(np.logical_and( + inds_ray_aux, inds_rng_aux = np.where( 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))) + 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') + 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_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] + azi = radar.azimuth["data"][ind_ray] + rng = radar.range["data"][ind_rng] return ind_ray, ind_rng, azi, rng @@ -1126,34 +1288,42 @@ def find_neighbour_gates(radar, azi, rng, delta_azi=None, delta_rng=None): """ # find gates close to lat lon point if delta_azi is None: - inds_ray = np.ma.arange(radar.azimuth['data'].size) + inds_ray = np.ma.arange(radar.azimuth["data"].size) else: azi_max = azi + delta_azi azi_min = azi - delta_azi - if azi_max > 360.: - azi_max -= 360. - if azi_min < 0.: - azi_min += 360. + 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] + 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] + 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) + 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] + 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., rmin=None, rmax=None): +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) @@ -1181,72 +1351,83 @@ def get_data_along_rng(radar, field_name, fix_elevations, fix_azimuths, """ if rmin is None: - rmin = 0. + rmin = 0.0 if rmax is None: - rmax = np.max(radar.range['data']) + rmax = np.max(radar.range["data"]) - rng_mask = np.logical_and( - radar.range['data'] >= rmin, radar.range['data'] <= rmax) + rng_mask = np.logical_and(radar.range["data"] >= rmin, radar.range["data"] <= rmax) - x = radar.range['data'][rng_mask] + x = radar.range["data"][rng_mask] xvals = [] yvals = [] valid_azi = [] valid_ele = [] - if radar.scan_type in ('ppi', 'vertical_pointing'): + if radar.scan_type in ("ppi", "vertical_pointing"): for ele, azi in zip(fix_elevations, fix_azimuths): - if radar.scan_type == 'vertical_pointing': + 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) + radar.fixed_angle["data"], ele, ang_tol=ang_tol + ) if ind_sweep is None: - warn('No elevation angle found for fix_elevation ' + str(ele)) + 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) + 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)) + 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]) + 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) + 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)) + 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) + 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)) + warn( + " No data found at azimuth " + + str(azi) + + " and elevation " + + str(ele) + ) continue - yvals.extend( - dataset_line.fields[field_name]['data'][0, rng_mask]) + 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]) + 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., ang_tol=1., azi_start=None, - azi_stop=None): +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) @@ -1276,26 +1457,24 @@ def get_data_along_azi(radar, field_name, fix_ranges, fix_elevations, """ if azi_start is None: - azi_start = np.min(radar.azimuth['data']) + azi_start = np.min(radar.azimuth["data"]) if azi_stop is None: - azi_stop = np.max(radar.azimuth['data']) + 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) + 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)) + 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 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)) + warn("No elevation angle found for fix_elevation " + str(ele)) continue new_dataset = radar.extract_sweeps([ind_sweep]) else: @@ -1303,29 +1482,37 @@ def get_data_along_azi(radar, field_name, fix_ranges, fix_elevations, 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)) + " 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) + 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]) + 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., ang_tol=1., ele_min=None, - ele_max=None): +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) @@ -1355,49 +1542,48 @@ def get_data_along_ele(radar, field_name, fix_ranges, fix_azimuths, """ if ele_min is None: - ele_min = np.min(radar.elevation['data']) + ele_min = np.min(radar.elevation["data"]) if ele_max is None: - ele_max = np.max(radar.elevation['data']) + 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) + 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)) + warn("No range gate found for fix_range " + str(rng)) continue - if radar.scan_type == 'ppi': + 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)) + " 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) + 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)) + 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]) + 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.): +def find_ang_index(ang_vec, ang, ang_tol=0.0): """ Find the angle index corresponding to a particular fixed angle @@ -1424,8 +1610,9 @@ def find_ang_index(ang_vec, ang, ang_tol=0.): return ind_ang -def _create_qvp_object(radar, field_names, qvp_type='qvp', start_time=None, - hmax=10000., hres=200.): +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 @@ -1458,45 +1645,45 @@ def _create_qvp_object(radar, field_names, qvp_type='qvp', start_time=None, 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') + qvp.fields[field_name]["data"] = np.array([], dtype="float64") # fixed radar objects parameters - qvp.range['data'] = np.arange(hmax / hres) * hres + hres / 2. - qvp.ngates = len(qvp.range['data']) + 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'] + qvp.time["units"] = radar.time["units"] else: - qvp.time['units'] = make_time_unit_str(start_time) + qvp.time["units"] = make_time_unit_str(start_time) - qvp.time['data'] = np.array([], dtype='float64') + qvp.time["data"] = np.array([], dtype="float64") qvp.scan_type = qvp_type - qvp.sweep_number['data'] = np.array([0], dtype='int32') + 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') + 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]]) + 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]]) + 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.], dtype='float64') + 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.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): +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 @@ -1528,44 +1715,42 @@ def _create_along_coord_object(radar, field_names, rng_values, fixed_angle, 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') + acoord.fields[field_name]["data"] = np.array([], dtype="float64") # fixed radar objects parameters - acoord.range['data'] = rng_values - acoord.ngates = len(acoord.range['data']) + acoord.range["data"] = rng_values + acoord.ngates = len(acoord.range["data"]) if start_time is None: - acoord.time['units'] = radar.time['units'] + acoord.time["units"] = radar.time["units"] else: - acoord.time['units'] = make_time_unit_str(start_time) + acoord.time["units"] = make_time_unit_str(start_time) - acoord.time['data'] = np.array([], dtype='float64') + acoord.time["data"] = np.array([], dtype="float64") acoord.scan_type = mode - acoord.sweep_number['data'] = np.array([0], dtype='int32') + 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') + 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]]) + 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.ray_angle_res["data"] = np.array([acoord.ray_angle_res["data"][0]]) - acoord.fixed_angle['data'] = np.array([fixed_angle], dtype='float64') + 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.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.): +def _update_qvp_metadata(qvp, ref_time, lon, lat, elev=90.0): """ updates a QVP object metadata with data from the current radar volume @@ -1582,23 +1767,22 @@ def _update_qvp_metadata(qvp, ref_time, lon, lat, elev=90.): 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 + 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. - qvp.elevation['data'] = ( - np.ones((qvp.nrays, ), dtype='float64') * elev) + 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)) + 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 @@ -1623,14 +1807,15 @@ def _update_along_coord_metadata(acoord, ref_time, elevation, azimuth): 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 + 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) + acoord.azimuth["data"] = np.append(acoord.azimuth["data"], azimuth) + acoord.elevation["data"] = np.append(acoord.elevation["data"], elevation) return acoord diff --git a/pyart/util/__init__.py b/pyart/util/__init__.py index 1a1674534de..b3591df6892 100644 --- a/pyart/util/__init__.py +++ b/pyart/util/__init__.py @@ -14,7 +14,7 @@ from .circular_stats import interval_std # noqa from .circular_stats import mean_of_two_angles # noqa from .circular_stats import mean_of_two_angles_deg # noqa -from .circular_stats import compute_directional_stats # noqa +from .circular_stats import compute_directional_stats # noqa from .columnsect import for_azimuth # noqa from .columnsect import get_column_rays # noqa from .columnsect import get_field_location # noqa @@ -33,7 +33,7 @@ determine_sweeps, subset_radar, to_vpt, - ma_broadcast_to + ma_broadcast_to, ) from .sigmath import ( # noqa angular_texture_2d, diff --git a/pyart/util/circular_stats.py b/pyart/util/circular_stats.py index b9a6003e1cb..c69359c978b 100644 --- a/pyart/util/circular_stats.py +++ b/pyart/util/circular_stats.py @@ -9,7 +9,8 @@ # https://en.wikipedia.org/wiki/Directional_statistics # https://en.wikipedia.org/wiki/Mean_of_circular_quantities -def compute_directional_stats(field, avg_type='mean', nvalid_min=1, axis=0): + +def compute_directional_stats(field, avg_type="mean", nvalid_min=1, axis=0): """ Computes the mean or the median along one of the axis (ray or range) @@ -31,7 +32,7 @@ def compute_directional_stats(field, avg_type='mean', nvalid_min=1, axis=0): nvalid : ndarray 1D The number of valid points used in the computation """ - if avg_type == 'mean': + if avg_type == "mean": values = np.ma.mean(field, axis=axis) else: values = np.ma.median(field, axis=axis) @@ -43,6 +44,7 @@ def compute_directional_stats(field, avg_type='mean', nvalid_min=1, axis=0): return values, nvalid + def mean_of_two_angles(angles1, angles2): """ Compute the element by element mean of two sets of angles. diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index 482068d022c..dd57eaad9f8 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -414,17 +414,17 @@ def subset_radar( if radar_aux.instrument_parameters is not None: if "nyquist_velocity" in radar_aux.instrument_parameters: - radar_aux.instrument_parameters["nyquist_velocity"][ - "data" - ] = radar_aux.instrument_parameters["nyquist_velocity"]["data"][ind_rays] + radar_aux.instrument_parameters["nyquist_velocity"]["data"] = ( + radar_aux.instrument_parameters["nyquist_velocity"]["data"][ind_rays] + ) if "pulse_width" in radar_aux.instrument_parameters: - radar_aux.instrument_parameters["pulse_width"][ - "data" - ] = radar_aux.instrument_parameters["pulse_width"]["data"][ind_rays] + radar_aux.instrument_parameters["pulse_width"]["data"] = ( + radar_aux.instrument_parameters["pulse_width"]["data"][ind_rays] + ) if "number_of_pulses" in radar_aux.instrument_parameters: - radar_aux.instrument_parameters["number_of_pulses"][ - "data" - ] = radar_aux.instrument_parameters["number_of_pulses"]["data"][ind_rays] + radar_aux.instrument_parameters["number_of_pulses"]["data"] = ( + radar_aux.instrument_parameters["number_of_pulses"]["data"][ind_rays] + ) # Get new fields if field_names is None: @@ -690,6 +690,7 @@ def image_mute_radar(radar, field, mute_field, mute_threshold, field_threshold=N radar.add_field("muted_" + field, muted_dict) return radar + def ma_broadcast_to(array, tup): """ Is used to guarantee that a masked array can be broadcasted without @@ -714,7 +715,7 @@ def ma_broadcast_to(array, tup): initial_fill_value = array.fill_value broadcasted_mask = np.broadcast_to(initial_mask, tup) return np.ma.array( - broadcasted_array, mask=broadcasted_mask, - fill_value=initial_fill_value) + broadcasted_array, mask=broadcasted_mask, fill_value=initial_fill_value + ) - return broadcasted_array \ No newline at end of file + return broadcasted_array diff --git a/pyart/xradar/accessor.py b/pyart/xradar/accessor.py index 8a374c88ac4..e79a1bcb150 100644 --- a/pyart/xradar/accessor.py +++ b/pyart/xradar/accessor.py @@ -3,7 +3,6 @@ """ - import copy import numpy as np diff --git a/tests/bridge/test_wradlib_bridge.py b/tests/bridge/test_wradlib_bridge.py index e777bdb2d74..d2d6a9cdda5 100644 --- a/tests/bridge/test_wradlib_bridge.py +++ b/tests/bridge/test_wradlib_bridge.py @@ -1,6 +1,5 @@ """ Unit Tests for Py-ART's io/mdv.py module. """ - import numpy as np import pytest diff --git a/tests/core/test_transforms.py b/tests/core/test_transforms.py index e07124f4337..d3a5ac20886 100644 --- a/tests/core/test_transforms.py +++ b/tests/core/test_transforms.py @@ -189,11 +189,14 @@ def test_cartesian_to_geographic_aeqd(): assert_almost_equal(lon, -100.0, 3) assert_almost_equal(lat, 40.0, 3) + def test_cartesian_to_antenna(): - r,az,el = transforms.cartesian_to_antenna(np.array([1000,3000,2000]), - np.array([1000,2000,3000]), - np.array([500,1000,1500])) + r, az, el = transforms.cartesian_to_antenna( + np.array([1000, 3000, 2000]), + np.array([1000, 2000, 3000]), + np.array([500, 1000, 1500]), + ) assert np.allclose(r, [1500, 3741.6573, 3905.1248]) - assert np.allclose(az, [45., 56.30993247, 33.69006753]) + assert np.allclose(az, [45.0, 56.30993247, 33.69006753]) assert np.allclose(el, [19.47122063, 15.50135957, 22.5885387]) diff --git a/tests/graph/test_cm.py b/tests/graph/test_cm.py index 562ef1a976b..4ad9841d5bd 100644 --- a/tests/graph/test_cm.py +++ b/tests/graph/test_cm.py @@ -1,6 +1,5 @@ """ Unit Tests for Py-ART's graph/cm.py module. """ - import matplotlib from pyart.graph import cm diff --git a/tests/graph/test_cm_colorblind.py b/tests/graph/test_cm_colorblind.py index 5060ff01433..e71ff0ccee5 100644 --- a/tests/graph/test_cm_colorblind.py +++ b/tests/graph/test_cm_colorblind.py @@ -1,6 +1,5 @@ """ Unit Tests for Py-ART's graph/cm_colorblind.py module. """ - import matplotlib from pyart.graph import cm_colorblind diff --git a/tests/graph/test_gridmapdisplay.py b/tests/graph/test_gridmapdisplay.py index b0ce20c2542..7068bb142cd 100644 --- a/tests/graph/test_gridmapdisplay.py +++ b/tests/graph/test_gridmapdisplay.py @@ -1,4 +1,5 @@ """ Unit Tests for Py-ART's graph/gridmapdisplay.py module. """ + # execute this script to create figure_gridmapdisplay_*.png files. # TODO use matplotlib's @image_comparison decorator to compare to file diff --git a/tests/graph/test_gridmapdisplay_basemap.py b/tests/graph/test_gridmapdisplay_basemap.py index 4e31d2095ab..34efc2e34de 100644 --- a/tests/graph/test_gridmapdisplay_basemap.py +++ b/tests/graph/test_gridmapdisplay_basemap.py @@ -1,4 +1,5 @@ """ Unit Tests for Py-ART's graph/gridmapdisplay_basemap.py module. """ + # execute this script to create figure_gridmapdisplay_*.png files. # TODO use matplotlib's @image_comparison decorator to compare to file diff --git a/tests/graph/test_radarmapdisplay.py b/tests/graph/test_radarmapdisplay.py index 9a5680b77e1..51d3976d78f 100644 --- a/tests/graph/test_radarmapdisplay.py +++ b/tests/graph/test_radarmapdisplay.py @@ -1,4 +1,5 @@ """ Unit Tests for Py-ART's graph/radarmapdisplay.py module. """ + # execute this script to create figure_plot_radar_display_*.png files. # TODO use matplotlib's @image_comparison decorator to compare to file diff --git a/tests/graph/test_radarmapdisplay_basemap.py b/tests/graph/test_radarmapdisplay_basemap.py index dcee1412c5b..a811010644d 100644 --- a/tests/graph/test_radarmapdisplay_basemap.py +++ b/tests/graph/test_radarmapdisplay_basemap.py @@ -1,4 +1,5 @@ """ Unit Tests for Py-ART's graph/radarmapdisplay.py module. """ + # execute this script to create figure_plot_radar_display_*.png files. # TODO use matplotlib's @image_comparison decorator to compare to file diff --git a/tests/retrieve/test_echo_class.py b/tests/retrieve/test_echo_class.py index fbe3ecf6f36..c1810455487 100644 --- a/tests/retrieve/test_echo_class.py +++ b/tests/retrieve/test_echo_class.py @@ -144,7 +144,7 @@ def test_hydroclass_semisupervised(): mass_centers = pyart.retrieve.echo_class._get_mass_centers(10e9) - hydro_nofreq = pyart.retrieve.hydroclass_semisupervised(radar_small)['hydro'] + hydro_nofreq = pyart.retrieve.hydroclass_semisupervised(radar_small)["hydro"] assert "units" in hydro_nofreq.keys() assert "standard_name" in hydro_nofreq.keys() @@ -156,17 +156,18 @@ def test_hydroclass_semisupervised(): assert_allclose(hydro_nofreq["data"][0][-5:], [3, 3, 3, 3, 3]) radar_small.instrument_parameters["frequency"] = {"data": np.array([5e9])} - hydro_freq = pyart.retrieve.hydroclass_semisupervised(radar_small)['hydro'] + hydro_freq = pyart.retrieve.hydroclass_semisupervised(radar_small)["hydro"] assert_allclose(hydro_freq["data"][0][0:5], [7, 7, 7, 7, 7]) assert_allclose(hydro_freq["data"][0][-5:], [3, 3, 3, 3, 3]) hydro = pyart.retrieve.hydroclass_semisupervised( radar_small, mass_centers=mass_centers - )['hydro'] + )["hydro"] assert_allclose(hydro["data"][0][0:5], [7, 7, 7, 7, 7]) assert_allclose(hydro["data"][0][-5:], [3, 3, 3, 3, 3]) + def test_data_limits_table(): dlimits_dict = pyart.retrieve.echo_class._data_limits_table() test_dict = { @@ -260,7 +261,7 @@ def test_assign_to_class(): ) mass_centers = pyart.retrieve.echo_class._get_mass_centers(10e9) - field_dict = {'Zh':zh, 'ZDR':zdr, 'KDP':kdp, 'RhoHV':rhohv, 'relH':relh} + field_dict = {"Zh": zh, "ZDR": zdr, "KDP": kdp, "RhoHV": rhohv, "relH": relh} hydroclass, _, _ = pyart.retrieve.echo_class._assign_to_class( field_dict, mass_centers diff --git a/tests/retrieve/test_qvp.py b/tests/retrieve/test_qvp.py index 211d62b0ac3..4cc4f832272 100644 --- a/tests/retrieve/test_qvp.py +++ b/tests/retrieve/test_qvp.py @@ -10,68 +10,175 @@ @pytest.fixture def test_radar(): test_radar = pyart.testing.make_empty_ppi_radar(1000, 360, 1) - test_radar.range['data'] *= 100 - test_radar.elevation['data'] = np.ones(test_radar.elevation['data'].shape) * 10. - test_radar.fixed_angle['data'] = np.array([10.]) - refl = 0.1 * np.arange(test_radar.ngates) + test_radar.range["data"] *= 100 + test_radar.elevation["data"] = np.ones(test_radar.elevation["data"].shape) * 10.0 + test_radar.fixed_angle["data"] = np.array([10.0]) + refl = 0.1 * np.arange(test_radar.ngates) refl = np.tile(refl, (test_radar.nrays, 1)) - test_radar.add_field('reflectivity', {'data': refl}) + test_radar.add_field("reflectivity", {"data": refl}) return test_radar + def test_compute_qvp(test_radar): - qvp = pyart.retrieve.compute_qvp(test_radar, ['reflectivity'], hmax = 10000) - qvp_refl = np.array([1.8999999999999875, 2.200000000000011, - 2.3999999999999853, 2.7000000000000175, 3.0, - 3.299999999999978, 3.599999999999994, - 3.9000000000000132, 4.200000000000027, 4.5, - 4.7000000000000295, 5.0, 5.299999999999969, - 5.599999999999963, 5.900000000000039, - 6.2000000000000135, 6.5, 6.7, 7.0, 7.300000000000015]) - - assert np.allclose(qvp.fields['reflectivity']['data'][0,10:30], qvp_refl) - assert len(qvp.range['data']) == 200 - assert np.all(qvp.range['data'] == np.arange(25, 10000, 50)) - assert qvp.azimuth['data'][0] == 0 + qvp = pyart.retrieve.compute_qvp(test_radar, ["reflectivity"], hmax=10000) + qvp_refl = np.array( + [ + 1.8999999999999875, + 2.200000000000011, + 2.3999999999999853, + 2.7000000000000175, + 3.0, + 3.299999999999978, + 3.599999999999994, + 3.9000000000000132, + 4.200000000000027, + 4.5, + 4.7000000000000295, + 5.0, + 5.299999999999969, + 5.599999999999963, + 5.900000000000039, + 6.2000000000000135, + 6.5, + 6.7, + 7.0, + 7.300000000000015, + ] + ) + + assert np.allclose(qvp.fields["reflectivity"]["data"][0, 10:30], qvp_refl) + assert len(qvp.range["data"]) == 200 + assert np.all(qvp.range["data"] == np.arange(25, 10000, 50)) + assert qvp.azimuth["data"][0] == 0 + def test_compute_rqvp(test_radar): - rqvp = pyart.retrieve.compute_rqvp(test_radar, ['reflectivity'], hres = 50, - hmax = 5000, rmax = 5000, weight_power = -1) - rqvp_refl = np.array([1.8999999999999875, 2.200000000000011, - 2.3999999999999853, 2.7000000000000175, 3.0, - 3.299999999999978, 3.599999999999994, - 3.9000000000000132, np.nan, np.nan, np.nan, np.nan, - np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, - np.nan]) - rqvp_refl_mask = np.array([False, False, False, False, False, False, False, False, - True, True, True, True, True, True, True, True, - True, True, True, True]) - - assert np.allclose(rqvp.fields['reflectivity']['data'][0,10:30], rqvp_refl) - assert np.allclose(rqvp.fields['reflectivity']['data'][0,10:30].mask, rqvp_refl_mask) - - assert len(rqvp.range['data']) == 100 - assert np.all(rqvp.range['data'] == np.arange(25, 5000, 50)) - assert rqvp.azimuth['data'][0] == 0 + rqvp = pyart.retrieve.compute_rqvp( + test_radar, ["reflectivity"], hres=50, hmax=5000, rmax=5000, weight_power=-1 + ) + rqvp_refl = np.array( + [ + 1.8999999999999875, + 2.200000000000011, + 2.3999999999999853, + 2.7000000000000175, + 3.0, + 3.299999999999978, + 3.599999999999994, + 3.9000000000000132, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + ] + ) + rqvp_refl_mask = np.array( + [ + False, + False, + False, + False, + False, + False, + False, + False, + True, + True, + True, + True, + True, + True, + True, + True, + True, + True, + True, + True, + ] + ) + + assert np.allclose(rqvp.fields["reflectivity"]["data"][0, 10:30], rqvp_refl) + assert np.allclose( + rqvp.fields["reflectivity"]["data"][0, 10:30].mask, rqvp_refl_mask + ) + + assert len(rqvp.range["data"]) == 100 + assert np.all(rqvp.range["data"] == np.arange(25, 5000, 50)) + assert rqvp.azimuth["data"][0] == 0 + def test_compute_evp(test_radar): - evp = pyart.retrieve.compute_evp(test_radar, ['reflectivity'], -97, 36, - delta_rng = 40000) + evp = pyart.retrieve.compute_evp( + test_radar, ["reflectivity"], -97, 36, delta_rng=40000 + ) + + evp_refl = np.array( + [ + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + 32.7, + 33.7, + 35.10000000000001, + 36.5, + 37.89999999999999, + 39.29999999999999, + 40.70000000000001, + ] + ) - evp_refl = np.array([np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, - np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, np.nan, - 32.7, 33.7, 35.10000000000001, 36.5, 37.89999999999999, - 39.29999999999999, 40.70000000000001]) + evp_refl_mask = np.array( + [ + [ + True, + True, + True, + True, + True, + True, + True, + True, + True, + True, + True, + True, + True, + False, + False, + False, + False, + False, + False, + False, + ] + ] + ) - evp_refl_mask = np.array([[ True, True, True, True, True, True, True, True, - True, True, True, True, True, False, False, False, - False, False, False, False]]) + assert np.allclose(evp.fields["reflectivity"]["data"][0, 10:30], evp_refl) + assert np.allclose(evp.fields["reflectivity"]["data"][0, 10:30].mask, evp_refl_mask) - assert np.allclose(evp.fields['reflectivity']['data'][0,10:30], evp_refl) - assert np.allclose(evp.fields['reflectivity']['data'][0,10:30].mask, evp_refl_mask) + assert len(evp.range["data"]) == 40 + assert np.all(evp.range["data"] == np.arange(125, 10000, 250)) + assert evp.azimuth["data"][0] == 0 - assert len(evp.range['data']) == 40 - assert np.all(evp.range['data'] == np.arange(125, 10000, 250)) - assert evp.azimuth['data'][0] == 0 def test_quasi_vertical_profile(): test_radar = pyart.testing.make_target_radar() diff --git a/tests/util/test_circular_stats.py b/tests/util/test_circular_stats.py index 6a8251f5f1a..e5574c719f8 100644 --- a/tests/util/test_circular_stats.py +++ b/tests/util/test_circular_stats.py @@ -7,33 +7,35 @@ def test_compute_directional_stats(): - field_1d = np.ma.array([1,1,3,4,5]) + field_1d = np.ma.array([1, 1, 3, 4, 5]) field_1d.mask = [True, False, False, False, False] field = np.tile(field_1d, (10, 1)) - mean, nvalid = circular_stats.compute_directional_stats(field, axis = 1) - median, nvalid = circular_stats.compute_directional_stats(field, axis = 1, - avg_type = 'median') + mean, nvalid = circular_stats.compute_directional_stats(field, axis=1) + median, nvalid = circular_stats.compute_directional_stats( + field, axis=1, avg_type="median" + ) assert np.all(mean == 3.25) assert np.all(median == 3.5) assert np.all(nvalid == 4) # Test with larger nb of nvalid_min - mean, nvalid = circular_stats.compute_directional_stats(field, axis = 1, - nvalid_min = 5) + mean, nvalid = circular_stats.compute_directional_stats(field, axis=1, nvalid_min=5) assert np.all(mean.mask) # Test other direction - mean, nvalid = circular_stats.compute_directional_stats(field, axis = 0) - median, nvalid = circular_stats.compute_directional_stats(field, axis = 0, - avg_type = 'median') + mean, nvalid = circular_stats.compute_directional_stats(field, axis=0) + median, nvalid = circular_stats.compute_directional_stats( + field, axis=0, avg_type="median" + ) - ref = np.ma.array([np.nan, 1, 3, 4, 5], mask = [1, 0 , 0, 0, 0]) + ref = np.ma.array([np.nan, 1, 3, 4, 5], mask=[1, 0, 0, 0, 0]) assert np.all(mean == ref) assert np.all(median == ref) assert np.all(nvalid == [0, 10, 10, 10, 10]) + def test_mean_of_two_angles(): mean = circular_stats.mean_of_two_angles(np.pi / 4 + 0.2, 7 * np.pi / 4) assert_almost_equal(mean, 0.10, 2) diff --git a/tests/util/test_radar_utils.py b/tests/util/test_radar_utils.py index 49d005eaf7a..0a7afbc61b5 100644 --- a/tests/util/test_radar_utils.py +++ b/tests/util/test_radar_utils.py @@ -104,15 +104,17 @@ def test_subset_radar(): assert radarcut.elevation["data"].max() <= ele_max assert list(radarcut.fields) == ["f1"] + def test_ma_broadcast_to(): buf = np.ma.zeros(5) - buf.mask = [1,1,0,0,0] + buf.mask = [1, 1, 0, 0, 0] buf_broad = pyart.util.radar_utils.ma_broadcast_to(buf, (10, 5)) - assert buf_broad.shape == (10,5) - assert buf_broad.mask.shape == (10,5) - expected_mask = np.tile(np.array([True, True, False, False, False]),[10,1]) + assert buf_broad.shape == (10, 5) + assert buf_broad.mask.shape == (10, 5) + expected_mask = np.tile(np.array([True, True, False, False, False]), [10, 1]) assert np.all(expected_mask == buf_broad.mask) + # read in example file radar = pyart.io.read_nexrad_archive(pyart.testing.NEXRAD_ARCHIVE_MSG31_FILE) @@ -137,4 +139,3 @@ def test_image_mute_radar(): n_nonmutez = np.sum(~mute_radar.fields["nonmuted_reflectivity"]["data"].mask) assert n_mutez + n_nonmutez == n_rhohv - From 3ff109b9686363926fca9f1c4423e957702f36de Mon Sep 17 00:00:00 2001 From: Daniel Wolfensberger Date: Wed, 28 Aug 2024 08:40:40 +0200 Subject: [PATCH 04/15] FIX: modified pyart examples to handle new hydroclass_semisupervised output --- examples/retrieve/plot_hydrometeor.py | 2 +- .../retrieve/plot_hydrometeor_class_x_band.py | 2 +- pyart/correct/phase_proc.py | 6 +++++- pyart/retrieve/qvp.py | 4 +--- pyart/util/radar_utils.py | 18 +++++++++--------- 5 files changed, 17 insertions(+), 15 deletions(-) diff --git a/examples/retrieve/plot_hydrometeor.py b/examples/retrieve/plot_hydrometeor.py index 43702ca962b..f3f062c0391 100644 --- a/examples/retrieve/plot_hydrometeor.py +++ b/examples/retrieve/plot_hydrometeor.py @@ -53,7 +53,7 @@ kdp_field="specific_differential_phase", rhv_field="uncorrected_cross_correlation_ratio", temp_field="temperature", -) +)["hydro"] radar.add_field("radar_echo_classification", hydro) diff --git a/examples/retrieve/plot_hydrometeor_class_x_band.py b/examples/retrieve/plot_hydrometeor_class_x_band.py index 1841562fd89..f168a18f7a5 100644 --- a/examples/retrieve/plot_hydrometeor_class_x_band.py +++ b/examples/retrieve/plot_hydrometeor_class_x_band.py @@ -98,7 +98,7 @@ rhv_field="RHOHV", temp_field="sounding_temperature", radar_freq=9.2e9, -) +)["hydro"] radar.add_field("hydro_classification", hydromet_class, replace_existing=True) diff --git a/pyart/correct/phase_proc.py b/pyart/correct/phase_proc.py index 2d76473035d..27b6901f1c9 100644 --- a/pyart/correct/phase_proc.py +++ b/pyart/correct/phase_proc.py @@ -108,7 +108,11 @@ def fzl_index(fzl, ranges, elevation, radar_height): p_r = 4.0 * Re / 3.0 z = ( radar_height - + (ranges**2 + p_r**2 + 2.0 * ranges * p_r * np.sin(elevation * np.pi / 180.0)) + + ( + ranges**2 + + p_r**2 + + 2.0 * ranges * p_r * np.sin(elevation * np.pi / 180.0) + ) ** 0.5 - p_r ) diff --git a/pyart/retrieve/qvp.py b/pyart/retrieve/qvp.py index 18ef9e52e12..50f9cc3034c 100644 --- a/pyart/retrieve/qvp.py +++ b/pyart/retrieve/qvp.py @@ -24,7 +24,7 @@ _update_qvp_metadata _update_along_coord_metadata -""" +""" "" "" "" "" "" "" "" "" "" "" "" "" "" "" from copy import deepcopy from warnings import warn @@ -106,7 +106,6 @@ def quasi_vertical_profile(radar, desired_angle=None, fields=None, gatefilter=No 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: @@ -431,7 +430,6 @@ def compute_rqvp( ) 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) diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index dd57eaad9f8..ba7063a98cb 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -414,17 +414,17 @@ def subset_radar( if radar_aux.instrument_parameters is not None: if "nyquist_velocity" in radar_aux.instrument_parameters: - radar_aux.instrument_parameters["nyquist_velocity"]["data"] = ( - radar_aux.instrument_parameters["nyquist_velocity"]["data"][ind_rays] - ) + radar_aux.instrument_parameters["nyquist_velocity"][ + "data" + ] = radar_aux.instrument_parameters["nyquist_velocity"]["data"][ind_rays] if "pulse_width" in radar_aux.instrument_parameters: - radar_aux.instrument_parameters["pulse_width"]["data"] = ( - radar_aux.instrument_parameters["pulse_width"]["data"][ind_rays] - ) + radar_aux.instrument_parameters["pulse_width"][ + "data" + ] = radar_aux.instrument_parameters["pulse_width"]["data"][ind_rays] if "number_of_pulses" in radar_aux.instrument_parameters: - radar_aux.instrument_parameters["number_of_pulses"]["data"] = ( - radar_aux.instrument_parameters["number_of_pulses"]["data"][ind_rays] - ) + radar_aux.instrument_parameters["number_of_pulses"][ + "data" + ] = radar_aux.instrument_parameters["number_of_pulses"]["data"][ind_rays] # Get new fields if field_names is None: From 6b24a9fb9b089b39ecf25b89e8c5f419ae66e65f Mon Sep 17 00:00:00 2001 From: Daniel Wolfensberger Date: Wed, 28 Aug 2024 09:36:12 +0200 Subject: [PATCH 05/15] FIX: fix in plot_hydrometeor_class_x_band.py example which was using wrong hydroclass labels --- .../retrieve/plot_hydrometeor_class_x_band.py | 25 ++++++++----------- 1 file changed, 11 insertions(+), 14 deletions(-) diff --git a/examples/retrieve/plot_hydrometeor_class_x_band.py b/examples/retrieve/plot_hydrometeor_class_x_band.py index f168a18f7a5..0018f7213d1 100644 --- a/examples/retrieve/plot_hydrometeor_class_x_band.py +++ b/examples/retrieve/plot_hydrometeor_class_x_band.py @@ -51,7 +51,7 @@ kdp_field="filtered_corrected_specific_diff_phase", rhv_field="RHOHV", temp_field="sounding_temperature", -) +)['hydro'] radar.add_field("hydro_classification", hydromet_class, replace_existing=True) @@ -69,7 +69,7 @@ rhv_field="RHOHV", temp_field="sounding_temperature", radar_freq=9.2e9, -) +)['hydro'] radar.add_field("hydro_classification", hydromet_class, replace_existing=True) @@ -97,8 +97,8 @@ kdp_field="filtered_corrected_specific_diff_phase", rhv_field="RHOHV", temp_field="sounding_temperature", - radar_freq=9.2e9, -)["hydro"] + radar_freq=9.2e9 +)['hydro'] radar.add_field("hydro_classification", hydromet_class, replace_existing=True) @@ -121,7 +121,6 @@ "Lime", "Yellow", "Red", - "Fuchsia", ] cmaphid = colors.ListedColormap(hid_colors) cmapmeth = colors.ListedColormap(hid_colors[0:6]) @@ -129,26 +128,24 @@ def adjust_fhc_colorbar_for_pyart(cb): - cb.set_ticks(np.arange(1.4, 10, 0.9)) + cb.set_ticks(np.arange(1.4, 9, 0.9)) cb.ax.set_yticklabels( [ - "Drizzle", - "Rain", - "Ice Crystals", "Aggregates", - "Wet Snow", + "Ice Crystals", + "Light rain", + "Rain", "Vertical Ice", - "LD Graupel", + "Wet snow", "HD Graupel", - "Hail", - "Big Drops", + "Melting hail", + "Dry hail or high-density graupel", ] ) cb.ax.set_ylabel("") cb.ax.tick_params(length=0) return cb - def adjust_meth_colorbar_for_pyart(cb, tropical=False): if not tropical: cb.set_ticks(np.arange(1.25, 5, 0.833)) From 7cb88e65ae3fb79f74d66ece4e329d70686511fe Mon Sep 17 00:00:00 2001 From: Daniel Wolfensberger Date: Wed, 28 Aug 2024 09:38:55 +0200 Subject: [PATCH 06/15] ENH: updated defaut_config.py to account for new outputs of hydroclass_semisupervised --- pyart/default_config.py | 80 +++++++++++++++++++++++++++++++++++++++++ 1 file changed, 80 insertions(+) diff --git a/pyart/default_config.py b/pyart/default_config.py index 876d654f475..dba449ce52a 100644 --- a/pyart/default_config.py +++ b/pyart/default_config.py @@ -91,6 +91,16 @@ rain_rate = "rain_rate" radar_estimated_rain_rate = "radar_estimated_rain_rate" radar_echo_classification = "radar_echo_classification" +hydroclass_entropy = "hydroclass_entropy" +proportion_AG = "proportion_AG" +proportion_CR = "proportion_CR" +proportion_LR = "proportion_LR" +proportion_RP = "proportion_RP" +proportion_RN = "proportion_RN" +proportion_VI = "proportion_VI" +proportion_WS = "proportion_WS" +proportion_MH = "proportion_MH" +proportion_IH = "proportion_IH" specific_attenuation = "specific_attenuation" specific_differential_attenuation = "specific_differential_attenuation" clutter_filter_power_removed = "clutter_filter_power_removed" @@ -640,6 +650,56 @@ "long_name": "Radar Echo classification", "coordinates": "elevation azimuth range", }, + hydroclass_entropy: { + "units": "-", + "standard_name": "hydroclass_entropy", + "long_name": "Semi-supervised hydrometeor classification entropy", + "coordinates": "elevation azimuth range"}, + proportion_AG: { + "units": "percent", + "standard_name": "proportion_AG", + "long_name": "Aggregates proportion", + "coordinates": "elevation azimuth range"}, + proportion_CR: { + "units": "percent", + "standard_name": "proportion_CR", + "long_name": "Crystals proportion", + "coordinates": "elevation azimuth range"}, + proportion_LR: { + "units": "percent", + "standard_name": "proportion_LR", + "long_name": "Light Rain proportion", + "coordinates": "elevation azimuth range"}, + proportion_RP: { + "units": "percent", + "standard_name": "proportion_RP", + "long_name": "Rimed particles proportion", + "coordinates": "elevation azimuth range"}, + proportion_RN: { + "units": "percent", + "standard_name": "proportion_RN", + "long_name": "Rain proportion", + "coordinates": "elevation azimuth range"}, + proportion_VI: { + "units": "percent", + "standard_name": "proportion_VI", + "long_name": "Vertical Ice Crystals proportion", + "coordinates": "elevation azimuth range"}, + proportion_WS: { + "units": "percent", + "standard_name": "proportion_WS", + "long_name": "Wet Snow proportion", + "coordinates": "elevation azimuth range"}, + proportion_MH: { + "units": "percent", + "standard_name": "proportion_MH", + "long_name": "Melting Hail proportion", + "coordinates": "elevation azimuth range"}, + proportion_IH: { + "units": "percent", + "standard_name": "proportion_IH", + "long_name": "Ice Hail proportion", + "coordinates": "elevation azimuth range"}, specific_attenuation: { "units": "dB/km", "standard_name": "specific_attenuation", @@ -1442,6 +1502,16 @@ def spectrum_width_limit(container=None, selection=0): rain_rate: "pyart_RRate11", radar_estimated_rain_rate: "pyart_RRate11", radar_echo_classification: "pyart_LangRainbow12", + hydroclass_entropy: "pyart_LangRainbow12", + proportion_AG: "pyart_LangRainbow12", + proportion_CR: "pyart_LangRainbow12", + proportion_LR: "pyart_LangRainbow12", + proportion_RP: "pyart_LangRainbow12", + proportion_RN: "pyart_LangRainbow12", + proportion_VI: "pyart_LangRainbow12", + proportion_WS: "pyart_LangRainbow12", + proportion_MH: "pyart_LangRainbow12", + proportion_IH: "pyart_LangRainbow12", specific_attenuation: "pyart_Carbone17", differential_phase_texture: "pyart_BlueBrown11", height: "pyart_SCook18", @@ -1492,6 +1562,16 @@ def spectrum_width_limit(container=None, selection=0): rain_rate: (0.0, 50.0), radar_estimated_rain_rate: (0.0, 50.0), radar_echo_classification: (0, 11), + hydroclass_entropy: (0., 1.), + proportion_AG: (0., 100.), + proportion_CR: (0., 100.), + proportion_LR: (0., 100.), + proportion_RP: (0., 100.), + proportion_RN: (0., 100.), + proportion_VI: (0., 100.), + proportion_WS: (0., 100.), + proportion_MH: (0., 100.), + proportion_IH: (0., 100.), specific_attenuation: (0.0, 10.0), differential_phase_texture: (0, 180.0), height: (0, 20000), From 50c60e2da72e2b7e95a7a7dd633ae198b99cf924 Mon Sep 17 00:00:00 2001 From: Daniel Wolfensberger Date: Wed, 28 Aug 2024 09:45:35 +0200 Subject: [PATCH 07/15] FIX: adapted hydroclass_semisupervised to account for specified radar_freq --- pyart/retrieve/echo_class.py | 31 ++++++++++++++++++++++--------- 1 file changed, 22 insertions(+), 9 deletions(-) diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py index 5cb61e73024..37f29213618 100644 --- a/pyart/retrieve/echo_class.py +++ b/pyart/retrieve/echo_class.py @@ -617,6 +617,7 @@ def hydroclass_semisupervised( iso0_field=None, hydro_field=None, entropy_field=None, + radar_freq=None, temp_ref="temperature", compute_entropy=False, output_distances=False, @@ -657,6 +658,13 @@ def hydroclass_semisupervised( Output. Field name which represents the hydrometeor class field. A value of None will use the default field name as defined in the Py-ART configuration file. + entropy_field : str + Output. Field name which represents the entropy class field. + A value of None will use the default field name as defined in the + Py-ART configuration file. + radar_freq : str, optional + Radar frequency in Hertz (Hz) used for classification. + This parameter will be ignored, if the radar object has frequency information. temp_ref : str the field use as reference for temperature. Can be either temperature or height_over_iso0 @@ -743,15 +751,20 @@ def hydroclass_semisupervised( # select the centroids as a function of frequency band if mass_centers is None: # assign coefficients according to radar frequency - if "frequency" in radar.instrument_parameters: - mass_centers = _get_mass_centers( - radar.instrument_parameters["frequency"]["data"][0] + if radar.instrument_parameters and "frequency" in radar.instrument_parameters: + frequency = radar.instrument_parameters["frequency"]["data"][0] + mass_centers = _get_mass_centers(frequency) + warn(f"Using radar frequency from instrument parameters: {frequency}") + elif radar_freq is not None: + mass_centers = _get_mass_centers(radar_freq) + warn( + f"Radar instrument parameters are empty. Using user-supplied radar frequency: {radar_freq}" ) else: mass_centers = _mass_centers_table()["C"] warn( - "Radar frequency unknown. " - + "Default coefficients for C band will be applied" + "Radar instrument parameters and radar_freq param are empty." + "So frequency is unknown. Default coefficients for C band will be applied." ) if hydro_field is None: @@ -983,8 +996,8 @@ def _assign_to_class( # Get hydrometeor class class_vec = dist.argsort(axis=0, fill_value=10e40) - hydroclass_ray = (class_vec[0, :] + 2).astype(np.uint8) - hydroclass_ray[mask] = 1 + hydroclass_ray = (class_vec[0, :] + 1).astype(np.uint8) + hydroclass_ray[mask] = 0 hydroclass[ray, :] = hydroclass_ray if t_vals is None: @@ -1089,8 +1102,8 @@ def _assign_to_class_scan( # Get hydrometeor class class_vec = dist.argsort(axis=-1, fill_value=10e40) - hydroclass = np.ma.asarray(class_vec[:, :, 0] + 2, dtype=np.uint8) - hydroclass[mask] = 1 + hydroclass = np.ma.asarray(class_vec[:, :, 0] + 1, dtype=np.uint8) + hydroclass[mask] = 0 if t_vals is not None: # Transform the distance using the coefficient of the dominant class From 8ed83c88c4abb9539093e8705a588fc5f49dd687 Mon Sep 17 00:00:00 2001 From: Daniel Wolfensberger Date: Wed, 28 Aug 2024 10:55:08 +0200 Subject: [PATCH 08/15] FIX: adapted hydroclass_semisupervised to account for specified radar_freq --- pyart/retrieve/echo_class.py | 2 +- tests/retrieve/test_qvp.py | 209 +++++++++++++++++++++++++++++++++++ 2 files changed, 210 insertions(+), 1 deletion(-) diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py index 37f29213618..63e201d9cd6 100644 --- a/pyart/retrieve/echo_class.py +++ b/pyart/retrieve/echo_class.py @@ -698,10 +698,10 @@ def hydroclass_semisupervised( - 9: Dry hail or high-density graupel if compute_entropy is True: - entropy : dict Shannon entropy of the hydrometeor demixing + if output_distances is True: propX : dict Proportion of a given hydrometeor class in the polarimetric decomposition of a radar volume diff --git a/tests/retrieve/test_qvp.py b/tests/retrieve/test_qvp.py index 4cc4f832272..2293ae56f7e 100644 --- a/tests/retrieve/test_qvp.py +++ b/tests/retrieve/test_qvp.py @@ -1,7 +1,10 @@ """ Unit Tests for Py-ART's retrieve/qvp.py module. """ +import datetime + import numpy as np import pytest +from netCDF4 import num2date from numpy.testing import assert_almost_equal import pyart @@ -353,3 +356,209 @@ def test_quasi_vertical_profile(): assert_almost_equal(qvp["height"], qvp_height, 3) assert_almost_equal(qvp["range"], qvp_range, 3) assert_almost_equal(qvp["reflectivity"], qvp_reflectivity, 3) + + +def test_project_to_vertical(): + z = np.array([100, 200, 500, 1000, 2000]) + data = np.ma.array(np.linspace(0, 50, len(z))) + znew = np.arange(0, 2300, 100) + data_out = pyart.retrieve.qvp.project_to_vertical(data, z, znew) + data_out2 = pyart.retrieve.qvp.project_to_vertical( + data, z, znew, interp_kind="nearest" + ) + ref_out = np.ma.array( + data=[ + np.nan, + 0.0, + 12.5, + np.nan, + np.nan, + 25.0, + np.nan, + np.nan, + np.nan, + np.nan, + 37.5, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + np.nan, + 50.0, + np.nan, + np.nan, + ], + mask=[ + True, + False, + False, + True, + True, + False, + True, + True, + True, + True, + False, + True, + True, + True, + True, + True, + True, + True, + True, + True, + False, + True, + True, + ], + ) + ref_out2 = np.ma.array( + data=[ + np.nan, + 0.0, + 12.5, + 12.5, + 25.0, + 25.0, + 25.0, + 25.0, + 37.5, + 37.5, + 37.5, + 37.5, + 37.5, + 37.5, + 37.5, + 37.5, + 50.0, + 50.0, + 50.0, + 50.0, + 50.0, + np.nan, + np.nan, + ], + mask=[ + True, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + False, + True, + True, + ], + ) + + assert np.all(data_out == ref_out) + assert np.all(data_out2 == ref_out2) + + +def test_find_rng_index(): + vec = np.arange(0, 1000, 10) + idx = pyart.retrieve.qvp.find_rng_index(vec, 140) + assert idx == 14 + + +def test_get_target_elevations(): + test_radar = pyart.testing.make_target_radar() + test_radar.elevation["data"][0] = 2 + target_elevations, el_tol = pyart.retrieve.qvp.get_target_elevations(test_radar) + assert target_elevations[-1] == 2 + assert np.all(target_elevations[0:-1] == 0.75) + assert el_tol == 0 + + +def test_find_nearest_gate(test_radar): + ind_ray, ind_rng, azi, rng = pyart.retrieve.qvp.find_nearest_gate( + test_radar, 36.4, -97.4 + ) + + assert ind_ray == 141.0 + assert ind_rng == 145.0 + assert azi == 141.0 + assert rng == 14514.514 + + +def test_find_neighbour_gates(test_radar): + inds_ray, inds_rng = pyart.retrieve.qvp.find_neighbour_gates( + test_radar, 141, 14514, delta_azi=3, delta_rng=200 + ) + + assert np.all(inds_ray == [139, 140, 141, 142, 143]) + assert np.all(inds_rng == [143, 144, 145, 146]) + + +def test_find_ang_index(): + ang_vec = np.arange(0, 360) + idx = pyart.retrieve.qvp.find_ang_index(ang_vec, 36) + assert idx == 36 + + +def test__create_qvp_object(test_radar): + qvp = pyart.retrieve.qvp._create_qvp_object( + test_radar, ["reflectivity"], hmax=5000, hres=500 + ) + assert np.all(qvp.range["data"] == np.arange(250, 5000, 500)) + assert "reflectivity" in qvp.fields + assert qvp.sweep_mode["data"] == ["qvp"] + assert qvp.fixed_angle["data"] == 10 + + +def test__create_along_coord_object(test_radar): + acoord = pyart.retrieve.qvp._create_along_coord_object( + test_radar, ["reflectivity"], np.arange(0, 10000, 100), 10, "ALONG_AZI" + ) + + assert np.all(acoord.range["data"] == np.arange(0, 10000, 100)) + assert acoord.sweep_mode["data"] == ["ALONG_AZI"] + + +def test__update_qvp_metadata(test_radar): + qvp = pyart.retrieve.qvp._create_qvp_object( + test_radar, ["reflectivity"], hmax=5000, hres=500 + ) + new_time = datetime.datetime(2024, 6, 10) + qvp = pyart.retrieve.qvp._update_qvp_metadata(qvp, new_time, 10, 45) + start_time = num2date(0, qvp.time["units"], qvp.time["calendar"]) + time_offset = (new_time - start_time).total_seconds() + + assert np.all(qvp.gate_longitude["data"] == 10) + assert np.all(qvp.gate_latitude["data"] == 45) + assert qvp.time["data"] == time_offset + + +def test__update_along_coord_metadata(test_radar): + acoord = pyart.retrieve.qvp._create_along_coord_object( + test_radar, ["reflectivity"], np.arange(0, 10000, 100), 10, "ALONG_AZI" + ) + new_time = datetime.datetime(2024, 6, 10) + acoord = pyart.retrieve.qvp._update_along_coord_metadata(acoord, new_time, 5, 200) + + start_time = num2date(0, acoord.time["units"], acoord.time["calendar"]) + time_offset = (new_time - start_time).total_seconds() + + assert acoord.time["data"] == time_offset + assert acoord.elevation["data"] == [5] + assert acoord.azimuth["data"] == [200] From cf58dc00efad7ec0b2b732a5ccfa9a421c47d1cf Mon Sep 17 00:00:00 2001 From: Daniel Wolfensberger Date: Wed, 28 Aug 2024 11:46:27 +0200 Subject: [PATCH 09/15] ADD: new tests for missing fcts --- tests/retrieve/test_echo_class.py | 54 +++++++++++++++++++++++-------- 1 file changed, 41 insertions(+), 13 deletions(-) diff --git a/tests/retrieve/test_echo_class.py b/tests/retrieve/test_echo_class.py index c1810455487..bc9469d3e09 100644 --- a/tests/retrieve/test_echo_class.py +++ b/tests/retrieve/test_echo_class.py @@ -144,29 +144,58 @@ def test_hydroclass_semisupervised(): mass_centers = pyart.retrieve.echo_class._get_mass_centers(10e9) - hydro_nofreq = pyart.retrieve.hydroclass_semisupervised(radar_small)["hydro"] - + hydro_nofreq = pyart.retrieve.hydroclass_semisupervised(radar_small)['hydro'] assert "units" in hydro_nofreq.keys() assert "standard_name" in hydro_nofreq.keys() assert "long_name" in hydro_nofreq.keys() assert "coordinates" in hydro_nofreq.keys() assert hydro_nofreq["data"].dtype == "uint8" - assert_allclose(hydro_nofreq["data"][0][0:5], [7, 7, 7, 7, 7]) - assert_allclose(hydro_nofreq["data"][0][-5:], [3, 3, 3, 3, 3]) + assert_allclose(hydro_nofreq["data"][0][0:5], [6, 6, 6, 6, 6]) + assert_allclose(hydro_nofreq["data"][0][-5:], [2, 2, 2, 2, 2]) radar_small.instrument_parameters["frequency"] = {"data": np.array([5e9])} - hydro_freq = pyart.retrieve.hydroclass_semisupervised(radar_small)["hydro"] - assert_allclose(hydro_freq["data"][0][0:5], [7, 7, 7, 7, 7]) - assert_allclose(hydro_freq["data"][0][-5:], [3, 3, 3, 3, 3]) + hydro_freq = pyart.retrieve.hydroclass_semisupervised(radar_small)['hydro'] + assert_allclose(hydro_freq["data"][0][0:5], [6, 6, 6, 6, 6]) + assert_allclose(hydro_freq["data"][0][-5:], [2, 2, 2, 2, 2]) hydro = pyart.retrieve.hydroclass_semisupervised( radar_small, mass_centers=mass_centers - )["hydro"] + )['hydro'] + assert_allclose(hydro["data"][0][0:5], [6, 6, 6, 6, 6]) + assert_allclose(hydro["data"][0][-5:], [2, 2, 2, 2, 2]) - assert_allclose(hydro["data"][0][0:5], [7, 7, 7, 7, 7]) - assert_allclose(hydro["data"][0][-5:], [3, 3, 3, 3, 3]) +def test_hydroclass_semisupervised_demixing(): + radar = pyart.io.read(pyart.testing.NEXRAD_ARCHIVE_MSG31_FILE) + temp_dict = pyart.config.get_metadata("temperature") + temp = np.empty(radar.fields["reflectivity"]["data"].shape) + temp[:] = -10.0 + temp_dict["data"] = temp + radar.add_field("temperature", temp_dict, replace_existing=True) + radar.fields["specific_differential_phase"] = radar.fields.pop("differential_phase") + radar_small = radar.extract_sweeps([0]) + + mass_centers = pyart.retrieve.echo_class._get_mass_centers(10e9) + + hydro = pyart.retrieve.hydroclass_semisupervised(radar_small, + mass_centers=mass_centers, + compute_entropy=True, + output_distances=True) + + assert_allclose(hydro['entropy']["data"][0][0:5].data, + [0.269866, 0.269866, 0.269866, 0.269866, 0.269866], + rtol = 1E-5) + assert_allclose(hydro['entropy']["data"][0][-5:].data, + [0.16527913, 0.16527913, 0.16527913, 0.165279136, 0.16527913], + rtol = 1E-5) + assert_allclose(hydro['proportion_CR']["data"][0][0:5].data, + [0.00713444, 0.00713444, 0.00713444, 0.00713444, 0.00713444], + rtol = 1E-5) + assert_allclose(hydro['proportion_CR']["data"][0][-5:].data, + [3.96782815e-07, 3.96782815e-07, 3.96782815e-07, 3.96782815e-07, + 3.96782815e-07], + rtol = 1E-5) def test_data_limits_table(): dlimits_dict = pyart.retrieve.echo_class._data_limits_table() @@ -261,13 +290,12 @@ def test_assign_to_class(): ) mass_centers = pyart.retrieve.echo_class._get_mass_centers(10e9) - field_dict = {"Zh": zh, "ZDR": zdr, "KDP": kdp, "RhoHV": rhohv, "relH": relh} + field_dict = {'Zh':zh, 'ZDR':zdr, 'KDP':kdp, 'RhoHV':rhohv, 'relH':relh} hydroclass, _, _ = pyart.retrieve.echo_class._assign_to_class( field_dict, mass_centers ) - assert_allclose(hydroclass[0], [8, 8, 8, 8, 8], atol=1e-7) - + assert_allclose(hydroclass[0], [7, 7, 7, 7, 7], atol=1e-7) def test_standardize(): kdp = np.array( From d75cf78a59b9d201453841008f54aff34da746c8 Mon Sep 17 00:00:00 2001 From: Daniel Wolfensberger Date: Wed, 28 Aug 2024 11:47:25 +0200 Subject: [PATCH 10/15] ENH: add support for more scan modes in qvp fcts --- pyart/retrieve/qvp.py | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/pyart/retrieve/qvp.py b/pyart/retrieve/qvp.py index 50f9cc3034c..18ef9e52e12 100644 --- a/pyart/retrieve/qvp.py +++ b/pyart/retrieve/qvp.py @@ -24,7 +24,7 @@ _update_qvp_metadata _update_along_coord_metadata -""" "" "" "" "" "" "" "" "" "" "" "" "" "" "" +""" from copy import deepcopy from warnings import warn @@ -106,6 +106,7 @@ def quasi_vertical_profile(radar, desired_angle=None, fields=None, gatefilter=No 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: @@ -430,6 +431,7 @@ def compute_rqvp( ) 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) From f6fe7e25b6ed7fc53a9107741f09d00fbe18379d Mon Sep 17 00:00:00 2001 From: Daniel Wolfensberger Date: Wed, 28 Aug 2024 11:57:58 +0200 Subject: [PATCH 11/15] FIX: docstring correction in cartesian_to_antenna --- pyart/core/transforms.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/pyart/core/transforms.py b/pyart/core/transforms.py index d37c2b3528d..4dff4b16e08 100644 --- a/pyart/core/transforms.py +++ b/pyart/core/transforms.py @@ -95,8 +95,6 @@ def cartesian_to_antenna(x, y, z): Azimuth angle of the radar in degrees. [-180., 180] elevations : array Elevation angle of the radar in degrees. - ke : float, optional - Effective radius scale factor """ ranges = np.sqrt(x**2.0 + y**2.0 + z**2.0) From 22b19e2c9bd778472aac7374fe25aac89c77c4c8 Mon Sep 17 00:00:00 2001 From: Daniel Wolfensberger Date: Wed, 28 Aug 2024 12:42:14 +0200 Subject: [PATCH 12/15] FIX: missing hydrometeor variables in default_config.py --- pyart/correct/phase_proc.py | 6 +--- pyart/default_config.py | 60 ++++++++++++++++++++++++------------- pyart/retrieve/qvp.py | 4 +-- pyart/util/radar_utils.py | 18 +++++------ 4 files changed, 51 insertions(+), 37 deletions(-) diff --git a/pyart/correct/phase_proc.py b/pyart/correct/phase_proc.py index 27b6901f1c9..2d76473035d 100644 --- a/pyart/correct/phase_proc.py +++ b/pyart/correct/phase_proc.py @@ -108,11 +108,7 @@ def fzl_index(fzl, ranges, elevation, radar_height): p_r = 4.0 * Re / 3.0 z = ( radar_height - + ( - ranges**2 - + p_r**2 - + 2.0 * ranges * p_r * np.sin(elevation * np.pi / 180.0) - ) + + (ranges**2 + p_r**2 + 2.0 * ranges * p_r * np.sin(elevation * np.pi / 180.0)) ** 0.5 - p_r ) diff --git a/pyart/default_config.py b/pyart/default_config.py index dba449ce52a..cba8690fc82 100644 --- a/pyart/default_config.py +++ b/pyart/default_config.py @@ -196,6 +196,16 @@ "rain_rate": rain_rate, "radar_estimated_rain_rate": radar_estimated_rain_rate, "radar_echo_classification": radar_echo_classification, + "hydroclass_entropy": hydroclass_entropy, + "proportion_AG": proportion_AG, + "proportion_CR": proportion_CR, + "proportion_LR": proportion_LR, + "proportion_RP": proportion_RP, + "proportion_RN": proportion_RN, + "proportion_VI": proportion_VI, + "proportion_WS": proportion_WS, + "proportion_MH": proportion_MH, + "proportion_IH": proportion_IH, "specific_attenuation": specific_attenuation, "differential_phase_texture": differential_phase_texture, "eastward_wind_component": eastward_wind_component, @@ -654,52 +664,62 @@ "units": "-", "standard_name": "hydroclass_entropy", "long_name": "Semi-supervised hydrometeor classification entropy", - "coordinates": "elevation azimuth range"}, + "coordinates": "elevation azimuth range", + }, proportion_AG: { "units": "percent", "standard_name": "proportion_AG", "long_name": "Aggregates proportion", - "coordinates": "elevation azimuth range"}, + "coordinates": "elevation azimuth range", + }, proportion_CR: { "units": "percent", "standard_name": "proportion_CR", "long_name": "Crystals proportion", - "coordinates": "elevation azimuth range"}, + "coordinates": "elevation azimuth range", + }, proportion_LR: { "units": "percent", "standard_name": "proportion_LR", "long_name": "Light Rain proportion", - "coordinates": "elevation azimuth range"}, + "coordinates": "elevation azimuth range", + }, proportion_RP: { "units": "percent", "standard_name": "proportion_RP", "long_name": "Rimed particles proportion", - "coordinates": "elevation azimuth range"}, + "coordinates": "elevation azimuth range", + }, proportion_RN: { "units": "percent", "standard_name": "proportion_RN", "long_name": "Rain proportion", - "coordinates": "elevation azimuth range"}, + "coordinates": "elevation azimuth range", + }, proportion_VI: { "units": "percent", "standard_name": "proportion_VI", "long_name": "Vertical Ice Crystals proportion", - "coordinates": "elevation azimuth range"}, + "coordinates": "elevation azimuth range", + }, proportion_WS: { "units": "percent", "standard_name": "proportion_WS", "long_name": "Wet Snow proportion", - "coordinates": "elevation azimuth range"}, + "coordinates": "elevation azimuth range", + }, proportion_MH: { "units": "percent", "standard_name": "proportion_MH", "long_name": "Melting Hail proportion", - "coordinates": "elevation azimuth range"}, + "coordinates": "elevation azimuth range", + }, proportion_IH: { "units": "percent", "standard_name": "proportion_IH", "long_name": "Ice Hail proportion", - "coordinates": "elevation azimuth range"}, + "coordinates": "elevation azimuth range", + }, specific_attenuation: { "units": "dB/km", "standard_name": "specific_attenuation", @@ -1562,16 +1582,16 @@ def spectrum_width_limit(container=None, selection=0): rain_rate: (0.0, 50.0), radar_estimated_rain_rate: (0.0, 50.0), radar_echo_classification: (0, 11), - hydroclass_entropy: (0., 1.), - proportion_AG: (0., 100.), - proportion_CR: (0., 100.), - proportion_LR: (0., 100.), - proportion_RP: (0., 100.), - proportion_RN: (0., 100.), - proportion_VI: (0., 100.), - proportion_WS: (0., 100.), - proportion_MH: (0., 100.), - proportion_IH: (0., 100.), + hydroclass_entropy: (0.0, 1.0), + proportion_AG: (0.0, 100.0), + proportion_CR: (0.0, 100.0), + proportion_LR: (0.0, 100.0), + proportion_RP: (0.0, 100.0), + proportion_RN: (0.0, 100.0), + proportion_VI: (0.0, 100.0), + proportion_WS: (0.0, 100.0), + proportion_MH: (0.0, 100.0), + proportion_IH: (0.0, 100.0), specific_attenuation: (0.0, 10.0), differential_phase_texture: (0, 180.0), height: (0, 20000), diff --git a/pyart/retrieve/qvp.py b/pyart/retrieve/qvp.py index 18ef9e52e12..11a621ad23b 100644 --- a/pyart/retrieve/qvp.py +++ b/pyart/retrieve/qvp.py @@ -24,7 +24,7 @@ _update_qvp_metadata _update_along_coord_metadata -""" +""" "" "" from copy import deepcopy from warnings import warn @@ -106,7 +106,6 @@ def quasi_vertical_profile(radar, desired_angle=None, fields=None, gatefilter=No 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: @@ -431,7 +430,6 @@ def compute_rqvp( ) 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) diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index ba7063a98cb..dd57eaad9f8 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -414,17 +414,17 @@ def subset_radar( if radar_aux.instrument_parameters is not None: if "nyquist_velocity" in radar_aux.instrument_parameters: - radar_aux.instrument_parameters["nyquist_velocity"][ - "data" - ] = radar_aux.instrument_parameters["nyquist_velocity"]["data"][ind_rays] + radar_aux.instrument_parameters["nyquist_velocity"]["data"] = ( + radar_aux.instrument_parameters["nyquist_velocity"]["data"][ind_rays] + ) if "pulse_width" in radar_aux.instrument_parameters: - radar_aux.instrument_parameters["pulse_width"][ - "data" - ] = radar_aux.instrument_parameters["pulse_width"]["data"][ind_rays] + radar_aux.instrument_parameters["pulse_width"]["data"] = ( + radar_aux.instrument_parameters["pulse_width"]["data"][ind_rays] + ) if "number_of_pulses" in radar_aux.instrument_parameters: - radar_aux.instrument_parameters["number_of_pulses"][ - "data" - ] = radar_aux.instrument_parameters["number_of_pulses"]["data"][ind_rays] + radar_aux.instrument_parameters["number_of_pulses"]["data"] = ( + radar_aux.instrument_parameters["number_of_pulses"]["data"][ind_rays] + ) # Get new fields if field_names is None: From fb22cea44e75dea4cf823c2bded3584460a1b892 Mon Sep 17 00:00:00 2001 From: Daniel Wolfensberger Date: Wed, 28 Aug 2024 13:01:53 +0200 Subject: [PATCH 13/15] fixes to make tests pass --- examples/plotting/radar-cross-section.ipynb | 5 +- examples/retrieve/column-example.ipynb | 7 +- .../retrieve/plot_hydrometeor_class_x_band.py | 9 ++- .../retrieve/wavelet_echo_class_example.ipynb | 10 +-- tests/retrieve/test_echo_class.py | 75 +++++++++---------- 5 files changed, 54 insertions(+), 52 deletions(-) diff --git a/examples/plotting/radar-cross-section.ipynb b/examples/plotting/radar-cross-section.ipynb index e5389e7d99b..4fe48a2ca0a 100644 --- a/examples/plotting/radar-cross-section.ipynb +++ b/examples/plotting/radar-cross-section.ipynb @@ -8,6 +8,7 @@ "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", + "\n", "import pyart\n", "from pyart.testing import get_test_data\n", "\n", @@ -344,7 +345,7 @@ "outputs": [], "source": [ "fig = plt.figure(figsize=(10, 8))\n", - "title = \"($^\\circ$N)\"\n", + "title = r\"($^\\circ$N)\"\n", "cross.corrected_reflectivity_horizontal.plot(\n", " cmap=\"pyart_HomeyerRainbow\",\n", " x=None,\n", @@ -362,7 +363,7 @@ "metadata": {}, "outputs": [], "source": [ - "from pyart.graph.common import generate_grid_time_begin, generate_field_name\n", + "from pyart.graph.common import generate_field_name, generate_grid_time_begin\n", "\n", "\n", "def generate_cross_section_title(grid, field, start, end):\n", diff --git a/examples/retrieve/column-example.ipynb b/examples/retrieve/column-example.ipynb index eeeca093759..5f6beb51fe1 100644 --- a/examples/retrieve/column-example.ipynb +++ b/examples/retrieve/column-example.ipynb @@ -16,10 +16,11 @@ "metadata": {}, "outputs": [], "source": [ - "import pyart\n", - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", "import cartopy.crs as ccrs\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "\n", + "import pyart\n", "from pyart.testing import get_test_data" ] }, diff --git a/examples/retrieve/plot_hydrometeor_class_x_band.py b/examples/retrieve/plot_hydrometeor_class_x_band.py index 0018f7213d1..2e8d40b7692 100644 --- a/examples/retrieve/plot_hydrometeor_class_x_band.py +++ b/examples/retrieve/plot_hydrometeor_class_x_band.py @@ -51,7 +51,7 @@ kdp_field="filtered_corrected_specific_diff_phase", rhv_field="RHOHV", temp_field="sounding_temperature", -)['hydro'] +)["hydro"] radar.add_field("hydro_classification", hydromet_class, replace_existing=True) @@ -69,7 +69,7 @@ rhv_field="RHOHV", temp_field="sounding_temperature", radar_freq=9.2e9, -)['hydro'] +)["hydro"] radar.add_field("hydro_classification", hydromet_class, replace_existing=True) @@ -97,8 +97,8 @@ kdp_field="filtered_corrected_specific_diff_phase", rhv_field="RHOHV", temp_field="sounding_temperature", - radar_freq=9.2e9 -)['hydro'] + radar_freq=9.2e9, +)["hydro"] radar.add_field("hydro_classification", hydromet_class, replace_existing=True) @@ -146,6 +146,7 @@ def adjust_fhc_colorbar_for_pyart(cb): cb.ax.tick_params(length=0) return cb + def adjust_meth_colorbar_for_pyart(cb, tropical=False): if not tropical: cb.set_ticks(np.arange(1.25, 5, 0.833)) diff --git a/examples/retrieve/wavelet_echo_class_example.ipynb b/examples/retrieve/wavelet_echo_class_example.ipynb index 1b92486415c..64c1870c0d8 100644 --- a/examples/retrieve/wavelet_echo_class_example.ipynb +++ b/examples/retrieve/wavelet_echo_class_example.ipynb @@ -74,8 +74,9 @@ } ], "source": [ - "import pyart\n", - "import numpy as np" + "import numpy as np\n", + "\n", + "import pyart" ] }, { @@ -281,10 +282,9 @@ "source": [ "# Required imports\n", "\n", - "import matplotlib.pyplot as plt\n", - "import matplotlib.colors as mcolors\n", "import cartopy.crs as ccrs\n", - "\n", + "import matplotlib.colors as mcolors\n", + "import matplotlib.pyplot as plt\n", "\n", "display = pyart.graph.GridMapDisplay(grid)\n", "\n", diff --git a/tests/retrieve/test_echo_class.py b/tests/retrieve/test_echo_class.py index bc9469d3e09..0c9ed1a2dee 100644 --- a/tests/retrieve/test_echo_class.py +++ b/tests/retrieve/test_echo_class.py @@ -144,7 +144,7 @@ def test_hydroclass_semisupervised(): mass_centers = pyart.retrieve.echo_class._get_mass_centers(10e9) - hydro_nofreq = pyart.retrieve.hydroclass_semisupervised(radar_small)['hydro'] + hydro_nofreq = pyart.retrieve.hydroclass_semisupervised(radar_small)["hydro"] assert "units" in hydro_nofreq.keys() assert "standard_name" in hydro_nofreq.keys() assert "long_name" in hydro_nofreq.keys() @@ -155,47 +155,45 @@ def test_hydroclass_semisupervised(): assert_allclose(hydro_nofreq["data"][0][-5:], [2, 2, 2, 2, 2]) radar_small.instrument_parameters["frequency"] = {"data": np.array([5e9])} - hydro_freq = pyart.retrieve.hydroclass_semisupervised(radar_small)['hydro'] + hydro_freq = pyart.retrieve.hydroclass_semisupervised(radar_small)["hydro"] assert_allclose(hydro_freq["data"][0][0:5], [6, 6, 6, 6, 6]) assert_allclose(hydro_freq["data"][0][-5:], [2, 2, 2, 2, 2]) hydro = pyart.retrieve.hydroclass_semisupervised( - radar_small, mass_centers=mass_centers - )['hydro'] - assert_allclose(hydro["data"][0][0:5], [6, 6, 6, 6, 6]) - assert_allclose(hydro["data"][0][-5:], [2, 2, 2, 2, 2]) - - -def test_hydroclass_semisupervised_demixing(): - radar = pyart.io.read(pyart.testing.NEXRAD_ARCHIVE_MSG31_FILE) - temp_dict = pyart.config.get_metadata("temperature") - temp = np.empty(radar.fields["reflectivity"]["data"].shape) - temp[:] = -10.0 - temp_dict["data"] = temp - radar.add_field("temperature", temp_dict, replace_existing=True) - radar.fields["specific_differential_phase"] = radar.fields.pop("differential_phase") - radar_small = radar.extract_sweeps([0]) - - mass_centers = pyart.retrieve.echo_class._get_mass_centers(10e9) + radar_small, + mass_centers=mass_centers, + compute_entropy=True, + output_distances=True, + ) + assert_allclose(hydro["hydro"]["data"][0][0:5], [6, 6, 6, 6, 6]) + assert_allclose(hydro["hydro"]["data"][0][-5:], [2, 2, 2, 2, 2]) + assert_allclose( + hydro["entropy"]["data"][0][0:5].data, + [0.269866, 0.269866, 0.269866, 0.269866, 0.269866], + rtol=1e-5, + ) + assert_allclose( + hydro["entropy"]["data"][0][-5:].data, + [0.16527913, 0.16527913, 0.16527913, 0.165279136, 0.16527913], + rtol=1e-5, + ) + assert_allclose( + hydro["proportion_CR"]["data"][0][0:5].data, + [0.00713444, 0.00713444, 0.00713444, 0.00713444, 0.00713444], + rtol=1e-5, + ) + assert_allclose( + hydro["proportion_CR"]["data"][0][-5:].data, + [ + 3.96782815e-07, + 3.96782815e-07, + 3.96782815e-07, + 3.96782815e-07, + 3.96782815e-07, + ], + rtol=1e-5, + ) - hydro = pyart.retrieve.hydroclass_semisupervised(radar_small, - mass_centers=mass_centers, - compute_entropy=True, - output_distances=True) - - assert_allclose(hydro['entropy']["data"][0][0:5].data, - [0.269866, 0.269866, 0.269866, 0.269866, 0.269866], - rtol = 1E-5) - assert_allclose(hydro['entropy']["data"][0][-5:].data, - [0.16527913, 0.16527913, 0.16527913, 0.165279136, 0.16527913], - rtol = 1E-5) - assert_allclose(hydro['proportion_CR']["data"][0][0:5].data, - [0.00713444, 0.00713444, 0.00713444, 0.00713444, 0.00713444], - rtol = 1E-5) - assert_allclose(hydro['proportion_CR']["data"][0][-5:].data, - [3.96782815e-07, 3.96782815e-07, 3.96782815e-07, 3.96782815e-07, - 3.96782815e-07], - rtol = 1E-5) def test_data_limits_table(): dlimits_dict = pyart.retrieve.echo_class._data_limits_table() @@ -290,13 +288,14 @@ def test_assign_to_class(): ) mass_centers = pyart.retrieve.echo_class._get_mass_centers(10e9) - field_dict = {'Zh':zh, 'ZDR':zdr, 'KDP':kdp, 'RhoHV':rhohv, 'relH':relh} + field_dict = {"Zh": zh, "ZDR": zdr, "KDP": kdp, "RhoHV": rhohv, "relH": relh} hydroclass, _, _ = pyart.retrieve.echo_class._assign_to_class( field_dict, mass_centers ) assert_allclose(hydroclass[0], [7, 7, 7, 7, 7], atol=1e-7) + def test_standardize(): kdp = np.array( ([180.5, 180.5, 180.5, 180.5, 180.5], [180.5, 180.5, 180.5, 180.5, 180.5]) From de64fba50139e7cef2b221c15f44fb2bb780c1b3 Mon Sep 17 00:00:00 2001 From: Daniel Wolfensberger Date: Wed, 28 Aug 2024 13:10:14 +0200 Subject: [PATCH 14/15] fixes to make tests pass --- pyart/correct/phase_proc.py | 6 +++++- pyart/retrieve/qvp.py | 2 +- pyart/util/radar_utils.py | 18 +++++++++--------- 3 files changed, 15 insertions(+), 11 deletions(-) diff --git a/pyart/correct/phase_proc.py b/pyart/correct/phase_proc.py index 2d76473035d..27b6901f1c9 100644 --- a/pyart/correct/phase_proc.py +++ b/pyart/correct/phase_proc.py @@ -108,7 +108,11 @@ def fzl_index(fzl, ranges, elevation, radar_height): p_r = 4.0 * Re / 3.0 z = ( radar_height - + (ranges**2 + p_r**2 + 2.0 * ranges * p_r * np.sin(elevation * np.pi / 180.0)) + + ( + ranges**2 + + p_r**2 + + 2.0 * ranges * p_r * np.sin(elevation * np.pi / 180.0) + ) ** 0.5 - p_r ) diff --git a/pyart/retrieve/qvp.py b/pyart/retrieve/qvp.py index 11a621ad23b..b4a6b4a63ee 100644 --- a/pyart/retrieve/qvp.py +++ b/pyart/retrieve/qvp.py @@ -24,7 +24,7 @@ _update_qvp_metadata _update_along_coord_metadata -""" "" "" +""" "" "" "" "" "" "" "" "" from copy import deepcopy from warnings import warn diff --git a/pyart/util/radar_utils.py b/pyart/util/radar_utils.py index dd57eaad9f8..ba7063a98cb 100644 --- a/pyart/util/radar_utils.py +++ b/pyart/util/radar_utils.py @@ -414,17 +414,17 @@ def subset_radar( if radar_aux.instrument_parameters is not None: if "nyquist_velocity" in radar_aux.instrument_parameters: - radar_aux.instrument_parameters["nyquist_velocity"]["data"] = ( - radar_aux.instrument_parameters["nyquist_velocity"]["data"][ind_rays] - ) + radar_aux.instrument_parameters["nyquist_velocity"][ + "data" + ] = radar_aux.instrument_parameters["nyquist_velocity"]["data"][ind_rays] if "pulse_width" in radar_aux.instrument_parameters: - radar_aux.instrument_parameters["pulse_width"]["data"] = ( - radar_aux.instrument_parameters["pulse_width"]["data"][ind_rays] - ) + radar_aux.instrument_parameters["pulse_width"][ + "data" + ] = radar_aux.instrument_parameters["pulse_width"]["data"][ind_rays] if "number_of_pulses" in radar_aux.instrument_parameters: - radar_aux.instrument_parameters["number_of_pulses"]["data"] = ( - radar_aux.instrument_parameters["number_of_pulses"]["data"][ind_rays] - ) + radar_aux.instrument_parameters["number_of_pulses"][ + "data" + ] = radar_aux.instrument_parameters["number_of_pulses"]["data"][ind_rays] # Get new fields if field_names is None: From 73ca678b564afec9addef32e6777ab85c76edc22 Mon Sep 17 00:00:00 2001 From: Daniel Wolfensberger Date: Wed, 28 Aug 2024 15:09:36 +0200 Subject: [PATCH 15/15] FIX: changed value in ZDR standardization from -5 to 5 to -1.5 to 5 according to Besic (2016) --- pyart/retrieve/echo_class.py | 13 ++++++++----- tests/retrieve/test_echo_class.py | 20 ++++++++++---------- 2 files changed, 18 insertions(+), 15 deletions(-) diff --git a/pyart/retrieve/echo_class.py b/pyart/retrieve/echo_class.py index 63e201d9cd6..327acb9c00b 100644 --- a/pyart/retrieve/echo_class.py +++ b/pyart/retrieve/echo_class.py @@ -848,11 +848,11 @@ def hydroclass_semisupervised( hydro.update({"_FillValue": 0}) labels = ["NC"] ticks = [1] - boundaries = [0.5, 1.5] + boundaries = [-0.5, 1.5] for i, hydro_name in enumerate(hydro_names): labels.append(hydro_name) - ticks.append(i + 2) - boundaries.append(i + 2.5) + ticks.append(i + 1) + boundaries.append(i + 1.5) hydro.update({"labels": labels, "ticks": ticks, "boundaries": boundaries}) fields_dict.update({"hydro": hydro}) @@ -910,7 +910,9 @@ def _standardize(data, field_name, mx=None, mn=None): data[data < -0.5] = -0.5 data = 10.0 * np.ma.log10(data + 0.6) elif field_name == "RhoHV": - data = 10.0 * np.ma.log10(1.0 - data) + # avoid infinite result + data[data > 1.0] = 1.0 + data = 10.0 * np.ma.log10(1.0000000000001 - data) mask = np.ma.getmaskarray(data) field_std = 2.0 * (data - mn) / (mx - mn) - 1.0 @@ -1228,9 +1230,10 @@ def _data_limits_table(): """ dlimits_dict = dict() dlimits_dict.update({"Zh": (60.0, -10.0)}) - dlimits_dict.update({"ZDR": (5.0, -5.0)}) + dlimits_dict.update({"ZDR": (5.0, -1.5)}) dlimits_dict.update({"KDP": (7.0, -10.0)}) dlimits_dict.update({"RhoHV": (-5.23, -50.0)}) + dlimits_dict.update({"RelH": (5000.0, -5000.0)}) return dlimits_dict diff --git a/tests/retrieve/test_echo_class.py b/tests/retrieve/test_echo_class.py index 0c9ed1a2dee..47348447477 100644 --- a/tests/retrieve/test_echo_class.py +++ b/tests/retrieve/test_echo_class.py @@ -169,27 +169,27 @@ def test_hydroclass_semisupervised(): assert_allclose(hydro["hydro"]["data"][0][-5:], [2, 2, 2, 2, 2]) assert_allclose( hydro["entropy"]["data"][0][0:5].data, - [0.269866, 0.269866, 0.269866, 0.269866, 0.269866], + [0.35945634, 0.35945634, 0.35945634, 0.35945634, 0.35945634], rtol=1e-5, ) assert_allclose( hydro["entropy"]["data"][0][-5:].data, - [0.16527913, 0.16527913, 0.16527913, 0.165279136, 0.16527913], + [0.232788, 0.232788, 0.232788, 0.232788, 0.232788], rtol=1e-5, ) assert_allclose( hydro["proportion_CR"]["data"][0][0:5].data, - [0.00713444, 0.00713444, 0.00713444, 0.00713444, 0.00713444], + [0.03524214, 0.03524214, 0.03524214, 0.03524214, 0.03524214], rtol=1e-5, ) assert_allclose( hydro["proportion_CR"]["data"][0][-5:].data, [ - 3.96782815e-07, - 3.96782815e-07, - 3.96782815e-07, - 3.96782815e-07, - 3.96782815e-07, + 8.336829723993246e-05, + 8.336829723993246e-05, + 8.336829723993246e-05, + 8.336829723993246e-05, + 8.336829723993246e-05, ], rtol=1e-5, ) @@ -199,11 +199,11 @@ def test_data_limits_table(): dlimits_dict = pyart.retrieve.echo_class._data_limits_table() test_dict = { "Zh": (60.0, -10.0), - "ZDR": (5.0, -5.0), + "ZDR": (5, -1.5), "KDP": (7.0, -10.0), "RhoHV": (-5.23, -50.0), + "RelH": (5000, -5000), } - assert isinstance(dlimits_dict, dict) assert dlimits_dict == test_dict