diff --git a/pyart/retrieve/comp_z.py b/pyart/retrieve/comp_z.py index 01634279143..95fbc9751b4 100644 --- a/pyart/retrieve/comp_z.py +++ b/pyart/retrieve/comp_z.py @@ -9,11 +9,12 @@ from netCDF4 import num2date from pandas import to_datetime from scipy.interpolate import interp2d +from scipy.interpolate import RectBivariateSpline from pyart.core import Radar -def composite_reflectivity(radar, field="reflectivity", gatefilter=None): +def composite_reflectivity(radar, field="reflectivity", gatefilter=None,same_nyquist=True,nyquist_vector_idx=0): """ Composite Reflectivity @@ -41,6 +42,14 @@ def composite_reflectivity(radar, field="reflectivity", gatefilter=None): gatefilter : GateFilter GateFilter instance. None will result in no gatefilter mask being applied to data. + same_nyquist: bool + During a volume scan (i.e., file) the PRF (nyqust velocity) can change. + This can create some odd artifacts at times with if data quality is low on certain scans. + To get around this, you can change the code to only take the max of scans with the same nyquist +/- 1 m/s. + Defult this will be on (True), but folks can turn this off. + nyquist_vector_idx: int + This integer works alongside the same_nyquist parameter. Do you want to match all the nyquists to sweep 0? use idx 0, sweep 1, use idx 1. + By default it is set to 0. Returns ------- @@ -58,6 +67,7 @@ def composite_reflectivity(radar, field="reflectivity", gatefilter=None): # Determine the lowest sweep (used for metadata and such) minimum_sweep = np.min(radar.sweep_number["data"]) + # loop over all measured sweeps for sweep in sorted(radar.sweep_number["data"]): # get start and stop index numbers @@ -66,6 +76,9 @@ def composite_reflectivity(radar, field="reflectivity", gatefilter=None): # grab radar data z = radar.get_field(sweep, field) z_dtype = z.dtype + + #get the nyquist, so we know which sweeps are the same sens. + nyquist = np.asarray([np.round(radar.get_nyquist_vel(sweep=sweep))]) # Use gatefilter if gatefilter is not None: @@ -102,18 +115,29 @@ def composite_reflectivity(radar, field="reflectivity", gatefilter=None): lat_0[-1, :] = lat_0[0, :] else: - # Configure the intperpolator - z_interpolator = interp2d(ranges, az, z, kind="linear") + # Configure the intperpolator + z_interpolator = RectBivariateSpline(az, ranges, z) + # Apply the interpolation - z = z_interpolator(ranges, azimuth_final) + z = z_interpolator(azimuth_final,ranges) # if first sweep, create new dim, otherwise concat them up if sweep == minimum_sweep: z_stack = copy.deepcopy(z[np.newaxis, :, :]) + nyquist_stack = copy.deepcopy(nyquist[np.newaxis,:]) else: z_stack = np.concatenate([z_stack, z[np.newaxis, :, :]]) + nyquist_stack = np.concatenate([nyquist_stack, nyquist[np.newaxis,:]]) + + #only stack up sweeps with the same nyquist + if same_nyquist: + left = np.where(nyquist_stack >= nyquist_stack[nyquist_vector_idx]-1)[0] + right = np.where(nyquist_stack <= nyquist_stack[nyquist_vector_idx]+1)[0] + same_ny = np.intersect1d(left,right) + z_stack = z_stack[same_ny] + # now that the stack is made, take max across vertical compz = z_stack.max(axis=0).astype(z_dtype) @@ -190,4 +214,4 @@ def composite_reflectivity(radar, field="reflectivity", gatefilter=None): azimuth, elevation, instrument_parameters=instrument_parameters, - ) + ) \ No newline at end of file