Region-Based Dealiasing is wrongly unfolding tornadic signatures

[SteepAtticStairs]

• Py-ART version: 1.12.2
• Python version: 3.10.9
• Operating System: MacOS Ventura

Description

I was trying to dealias velocity data from the 2011 Super Outbreak. I was using this radar file (links to a direct download from the AWS S3 Bucket):
KBMX20110427_221945_V03.gz

What I Did

I used this python code to dealias and plot the radar file:

import matplotlib.pyplot as plt
import pyart
from pyart.correct import dealias_region_based


# read in file
radar = pyart.io.read_nexrad_archive('KBMX20110427_221945_V03.gz')

radar = radar.extract_sweeps([1]) # extract velocity sweep

theCorrVel = dealias_region_based(radar)
radar.add_field("dealiased_velocity", theCorrVel, replace_existing=True)

rmin = radar.fields['dealiased_velocity']['valid_min']
rmax = radar.fields['dealiased_velocity']['valid_max']

fig = plt.figure(figsize=[10, 7])
display = pyart.graph.RadarDisplay(radar)
display.plot('dealiased_velocity', vmin=rmin, vmax=rmax, cmap='pyart_balance')
plt.show()

However, this incorrectly unfolds the velocity data. Below are some screenshots detailing this:

This is clearly wrong, as the actual dealiasing should have two sections with massive gate-to-gate shear. However, the aliased data in that couplet do not have clear regions of folding as is usual, so that is giving the program some issues.

[SteepAtticStairs]

Another radar file from the same outbreak produces a similar error (again, this links to a direct download):
KBMX20110427_213353_V03.gz

[mgrover1]

@SteepAtticStairs thank you for expressing your concern with the algorithm. Dealiasing velocities around velocity couplets is a historically difficult task, and something that the region-based dealiaser does indeed struggle with. I encourage you check out this discussion thread:
https://openradar.discourse.group/t/over-aggressive-filtration/160

And experiment with with modifying the skip_between_rays and skip_along_ray parameters.

Often times in these regions, manual dealiasing is required using tools such as Hawkeye, which is a tool in the LROSE software ecosystem http://lrose.net/howtorun_HawkEye.html .

[SteepAtticStairs]

@mgrover1 Thanks for the response.

I have indeed experimented with those parameters, but they only seem to help when aliased regions beyond a range folded circle are not identified. For TVSs, I believe the issue is with the nyquist velocity.

This comment will concern the latter radar file that I brought up. In the aliased TVS that I focused on, the two regions of the shear are not clearly separated, even within the nyquist interval. This suggests to me that when the code is splitting up the nyquist interval to identify regions, the regions of the TVS are not different enough to be processed in a different chunk of the nyquist interval.

This lead me to experiment with the interval_splits parameter, but that too didn't fix the issue. Would splitting the nyquist interval into non-linear sections provide any meaningful difference? e.g. not splitting it into equally sized chunks.

ref: pyart/correct/region_dealias.py#L306

I'm almost certain that the dividing of the nyquist interval is the issue here, because other algorithms such as 4DD correctly dealias the TVS - but I know pyart is aiming to migrate away from the TRMM RSL and 4DD.

I hope this has made sense.

[mgrover1]

Yeah, there are pros/cons to each of these approaches, where the region based works generally throughout most of the domain, but struggles in these smaller regions with TVS features. I encourage you to keep experimenting, and adding an example showing how to tune this would be helpful for the documentation.

[mgrover1]

@SteepAtticStairs - coming back to this, I played around with the interval_splits parameter a bit, which did improve it slightly, but there is still an issue with the outbounds.

import matplotlib.pyplot as plt
import pyart
from pyart.correct import dealias_region_based


# read in file
radar = pyart.io.read_nexrad_archive('KBMX20110427_221945_V03.gz')

radar = radar.extract_sweeps([1]) # extract velocity sweep

theCorrVel = dealias_region_based(radar,
                                  interval_splits=30,
                                  centered=True)
radar.add_field("dealiased_velocity", theCorrVel, replace_existing=True)

rmin = radar.fields['dealiased_velocity']['valid_min']
rmax = radar.fields['dealiased_velocity']['valid_max']

display = pyart.graph.RadarDisplay(radar)
display.plot('dealiased_velocity', vmin=rmin, vmax=rmax, cmap='twilight_shifted')

plt.xlim(-70, -50)
plt.ylim(0, 20)
plt.ylim()
plt.show()

[SteepAtticStairs]

@mgrover1 Thanks for continuing to look into this.

It's interesting that simply adding more sections of the nyquist interval to region search has a slight difference for the better. However, I assume that scipy.ndimage.label isn't the quickest routine in the world, but the performance is still quite good for 30 executions.

I've been experimenting with non-linear nyquist interval splits, to some slight success, but it isn't anything that pyart could come up with on the fly - it seems specific to each dataset.

While I continue looking at this, I'll draw your attention to two new radar files with TVSs - one that doesn't dealias correctly, and one that does.

KSGF20110522_224838_V03.gz
This is from the 2011 Joplin tornado - this image is after dealiasing.

plt.xlim(-105, -85)
plt.ylim(-27, -8)

KTLX20130520_201643_V06.gz
This is from the 2013 Moore tornado. Surprisingly, a rather violent velocity couplet is dealiased pretty much perfectly.

plt.xlim(-33, -13)
plt.ylim(-11, 8)

[scollis]

Almost sounds like the beginnings of a study and a paper! This is all very interesting and I think it would be great to clone this conversation over https://openradar.discourse.group/ so it can get some wider attention (not saying this is not the appropriate place for discussing it here)

[scollis]

@rcjackson Could be good to test with code from @jmkurdzo
https://github.com/mit-ll/unet-vda

[mgrover1]

This is really interesting @SteepAtticStairs - would you be willing to work together on a blog post exploring this in detail, digging into what worked and didn't work? We can format it as a notebook! I think this will be helpful for the community.

[SteepAtticStairs]

@mgrover1 @scollis

I indeed would be interested in working on a blog post! You'd have to tell me the specifics on how that would be done, because I haven't really done something like this before. Also, I'm not sure how much help I'd be, since I know that a lot of people working on the pyart project have been in the meteorology field for a while, and have degrees and doctorates and all of that good stuff - which cannot be said about me.

It would also be helpful to clone this to Open Radar Discourse, but again, I'm not sure how that would happen, so I'd appreciate some help in that regard.

However, while we're doing this, I've discovered some parameters that successfully dealiased one of the problematic radar files mentioned before. I did this by running dealias_region_based twice, with a different interval_splits parameter for the second run. Here's the result:

import matplotlib.pyplot as plt
import pyart
import numpy as np

radar = pyart.io.read_nexrad_archive('KBMX20110427_213353_V03.gz')
radar = radar.extract_sweeps([1]) # extract velocity sweep

corrected_vel = pyart.correct.dealias_region_based(radar)
radar.add_field('velocity', corrected_vel, replace_existing=True)
corrected_vel = pyart.correct.dealias_region_based(radar, interval_splits=11)
radar.add_field('dealiased_velocity', corrected_vel, replace_existing=True)

rmin = radar.fields['dealiased_velocity']['valid_min']
rmax = radar.fields['dealiased_velocity']['valid_max']

fig = plt.figure(figsize=[10, 7])
display = pyart.graph.RadarDisplay(radar)
display.plot('dealiased_velocity', vmin=rmin, vmax=rmax, cmap='pyart_balance')

plt.xlim(-73, -53)
plt.ylim(42, 62)
plt.ylim()
plt.show()

[mgrover1]

@SteepAtticStairs this is great! Let's start by working on a notebook - do you use jupyter notebooks usually to develop these tests?

[SteepAtticStairs]

@mgrover1 Unfortunately not, I simply write it all down in a main.py file. However, I would be able to translate it to a notebook easily if needed.

[SteepAtticStairs]

@mgrover1 I've made a little progress, but I'm not sure how useful this will be.

My approach in this case was to:

1. Calculate the velocity texture of the entire field
2. Using that texture, filter out tornadic signatures
3. Region dealias the field without the tornadoes
4. Using that same texture (and the reflectivity field), isolate the tornadic signatures
5. Region dealias the isolated tornadoes
6. Step 5 will create artifacts, since dealias_region_based isn't great for TVSs, as has been established, so correct any errors by running dealias_unwrap_phase on the already region dealiased isolated tornadoes
7. Combine the tornadic and non-tornadic fields

This worked well for some radar files, and not so well for others. Maybe others can play around with this and fine-tune it? Or perhaps there is an easier way that I'm not thinking of.

import matplotlib.pyplot as plt
import pyart
import numpy as np

# read in file
radar = pyart.io.read_nexrad_archive('KBMX20110427_221945_V03.gz')
radar = radar.extract_sweeps([1]) # extract velocity sweep

# Calculate the velocity texture
vel_texture = pyart.retrieve.calculate_velocity_texture(radar, vel_field='velocity')
radar.add_field('vel_texture', vel_texture, replace_existing=True)

# Create a GateFilter to filter out tornadic signatures
gfilter_nontornadic = pyart.filters.GateFilter(radar)
# Tornadoes generally have low velocity texture values
gfilter_nontornadic.exclude_above('vel_texture', 5)

# Dealias the velocity field, filtering out all tornadic signatures
corrected_vel = pyart.correct.dealias_region_based(radar, gatefilter=gfilter_nontornadic)
radar.add_field('dealiased_nontornadic', corrected_vel, replace_existing=True)


# Create a GateFilter isolating tornadic signatures
gfilter_tornadic = pyart.filters.GateFilter(radar)
# Tornadoes AND noise have low velocity textures,
# but only tornadoes have a high reflectivity value,
# which is what we're trying to pinpoint here
gfilter_tornadic.exclude_below('reflectivity', 30)
gfilter_tornadic.exclude_below('vel_texture', 5)

# Region dealias the isolated tornadic signatures
corrected_vel = pyart.correct.dealias_region_based(radar)
radar.add_field('temp_dealiased_velocity', corrected_vel, replace_existing=True)
# Correct any artifacts in the isolated dealiased tornadoes by using "dealias_unwrap_phase"
corrected_vel = pyart.correct.dealias_unwrap_phase(radar, vel_field='temp_dealiased_velocity', gatefilter=gfilter_tornadic)
radar.add_field('dealiased_tornadic', corrected_vel, replace_existing=True)


# Combine the dealiased nontornadic and tornadic velocity fields
nontornadic = radar.fields['dealiased_nontornadic']['data']
tornadic = radar.fields['dealiased_tornadic']['data']

dealiased = np.ma.masked_invalid(nontornadic)
dealiased = np.ma.where(dealiased.mask, tornadic, dealiased)
radar.add_field_like('dealiased_nontornadic', 'dealiased', dealiased)


# Plot
field_to_plot = 'dealiased'

rmin = radar.fields[field_to_plot]['valid_min']
rmax = radar.fields[field_to_plot]['valid_max']

fig = plt.figure(figsize=[10, 7])
display = pyart.graph.RadarDisplay(radar)
display.plot(field_to_plot, vmin=rmin, vmax=rmax, cmap='pyart_balance')

plt.xlim(-73, -53)
plt.ylim(-3, 17)
plt.show()

[scollis]

That is a sweet TVS

[mgrover1]

This is really neat - I can draft up the notebook, and ping you once it is ready to look and at and you can suggest edits!!