Skip to content

Commit

Permalink
changes to comp_z
Browse files Browse the repository at this point in the history
  • Loading branch information
dopplerchase committed Apr 15, 2024
1 parent 4d694e3 commit ede33b0
Showing 1 changed file with 29 additions and 5 deletions.
34 changes: 29 additions & 5 deletions pyart/retrieve/comp_z.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
-------
Expand All @@ -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
Expand All @@ -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:
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -190,4 +214,4 @@ def composite_reflectivity(radar, field="reflectivity", gatefilter=None):
azimuth,
elevation,
instrument_parameters=instrument_parameters,
)
)

0 comments on commit ede33b0

Please sign in to comment.