-
Notifications
You must be signed in to change notification settings - Fork 265
New issue
Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.
By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.
Already on GitHub? Sign in to your account
Velocity unfolding #119
Comments
@kirknorth I'm not seeing these issues if I dealias using the original fields in the file. Can you share your processing script? Are you using the noise corrected reflectivity as the I want to have a good look at the FourDD algorithm in the future to implement some improvements. If I had some good example of problems to fix at the same time that would be great! I don't have the time that this requires right now, but hopefully after the AMS meeting in February. Here are my results, the interpolatedsonde file I'm using has a history attribute of: import matplotlib
matplotlib.rcParams['axes.titlesize'] = 'small'
import matplotlib.pyplot as plt
import netCDF4
import pyart
SOND_NAME = 'sgpinterpolatedsondeC1.c1.20110510.000000.cdf'
RADAR_NAME = 'sur/20110520/105235.mdv'
# read in the data
radar = pyart.io.read_mdv(RADAR_NAME)
# find and extract sonde data
target = netCDF4.num2date(radar.time['data'][0], radar.time['units'])
interp_sounde = netCDF4.Dataset(SOND_NAME)
t = pyart.correct.find_time_in_interp_sonde(interp_sounde, target)
height, speed, direction = t
# perform dealiasing
dealias_data = pyart.correct.dealias_fourdd(radar, height * 1000.0, speed,
direction, target)
radar.add_field('corrected_velocity', dealias_data)
# create a plots of velocity
fig = plt.figure(figsize=(12, 15))
display = pyart.graph.RadarDisplay(radar)
for i, sweep_i in enumerate([0, 1, 3, 6, 10, 14]):
ax = fig.add_subplot(3,2,i+1)
display.plot('velocity', sweep_i, vmin=-20, vmax=20, ax=ax,
colorbar_label='m/s', fig=fig)
display.set_limits((-120, 120), (-120, 120))
plt.tight_layout()
fig.savefig('velocity.png')
# plot corrected velocity
fig = plt.figure(figsize=(12, 15))
display = pyart.graph.RadarDisplay(radar)
for i, sweep_i in enumerate([0, 1, 3, 6, 10, 14]):
ax = fig.add_subplot(3,2,i+1)
display.plot('corrected_velocity', sweep_i, vmin=-20, vmax=20, ax=ax,
colorbar_label='m/s')
display.set_limits((-120, 120), (-120, 120))
plt.tight_layout()
fig.savefig('corrected_velocity.png') |
@jjhelmus, let me reprocess this file again. However, the things different about my processing and yours are the following:
And no, I did not use the noise corrected reflectivity field during dealiasing. I know this because I do that correction after calling |
@kirknorth I'm think this shows that something is wrong with the initial guess that Py-ART is providing to the FourDD algorithm, or that the algorithms is very unstable with regards to the initial guess. I'll look into it when I have a time to really dig into the code. Thanks for pointing this out, it helps to have a "bad" case when trying to make improvements. Just for completeness here is the new script: import matplotlib
matplotlib.rcParams['axes.titlesize'] = 'small'
import matplotlib.pyplot as plt
import netCDF4
import pyart
SOND_NAME = 'sgpinterpolatedsondeC1.c1.20110520.000000.cdf'
RADAR_NAME = 'sur/20110520/105235.mdv'
# read in the data
radar = pyart.io.read_mdv(RADAR_NAME)
# find and extract sonde data
target = netCDF4.num2date(radar.time['data'][0], radar.time['units'])
interp_sounde = netCDF4.Dataset(SOND_NAME)
t = pyart.correct.find_time_in_interp_sonde(interp_sounde, target)
height, speed, direction = t
# perform dealiasing
dealias_data = pyart.correct.dealias_fourdd(radar, height * 1000.0, speed,
direction, target,
keep_original=True)
radar.add_field('corrected_velocity', dealias_data)
# create a plots of velocity
fig = plt.figure(figsize=(12, 15))
display = pyart.graph.RadarDisplay(radar)
for i, sweep_i in enumerate([0, 1, 3, 6, 10, 14]):
ax = fig.add_subplot(3,2,i+1)
display.plot('velocity', sweep_i, vmin=-20, vmax=20, ax=ax,
colorbar_label='m/s', fig=fig)
display.set_limits((-120, 120), (-120, 120))
plt.tight_layout()
fig.savefig('velocity.png')
# plot corrected velocity
fig = plt.figure(figsize=(12, 15))
display = pyart.graph.RadarDisplay(radar)
for i, sweep_i in enumerate([0, 1, 3, 6, 10, 14]):
ax = fig.add_subplot(3,2,i+1)
display.plot('corrected_velocity', sweep_i, vmin=-20, vmax=20, ax=ax,
colorbar_label='m/s', mask_outside=False)
display.set_limits((-120, 120), (-120, 120))
plt.tight_layout()
fig.savefig('corrected_velocity.png') And if we expand the upper corrected velocity limit to 60 m/s: display.plot('corrected_velocity', sweep_i, vmin=-20, vmax=60, ax=ax, We can see that badly unfolded regions seem to be being assigned velocities in the 20-60 m/s range: |
For what is is worth I get good dealiasing if I set the speed and direction to zero for all points except the first height used by FourDD: ...
height, speed, direction = t
speed[2:284] = 0
direction[2:284] = 0
... I think there is something funny going on with the initial wind parameters from the sonde and how Py-ART or FourDD works with them. |
Yes.. Issue may be with the sonde -sent from a mobile device-
|
@jjhelmus we're certainly not short on examples where the velocity unfolding behaves oddly when using the ARM merged sounding product as an initial guess. Just for reference, here's another example, this time also showing how these artifacts end up in the gridded product. |
I meant to bring this up on the call today.. I actually think there may me issues with merged sounding.. Scott Collis On Jan 27, 2014, at 7:01 PM, Kirk North notifications@github.com wrote:
|
@scollis there very well could be issues with the merged sounding product. I can probably take a look at that tomorrow. However, I also checked the original (and I mean the original!) MMCG files you gave me way back when, and they show similar artifacts. And I don't believe you were using the merged sounding product as the initial guess back then... |
Yes.. Depends how original we speak.. Latest do use Merged Sonde from ARM.. As a community we need to investigate this.. Scott Collis On Jan 27, 2014, at 7:08 PM, Kirk North notifications@github.com wrote:
|
@kirknorth Just FYI, it would probably be good to expand your vmin/vmax limits when plotting the correct velocities. The original Doppler limits are +/- 20 m/s so any folded signals would be outside this range. Since dealiasing tried to unfold these velocities outside these limits are expected. +/- 40 or +/- 60 seem to work well when there is a single or double fold which seems to be the most common case. |
@jjhelmus, in the first couple of lines of refl_array = radar.fields[refl_field]['data'].reshape(nshape).astype('float32')
refl_array[np.where(refl_array == fill_value)] = rsl_badval
refl_array[np.where(np.isnan(refl_array))] = rsl_badval
refl_volume = _rsl_interface.create_volume(refl_array, 0) The Would the fix be to simply use the |
Okay so with the new (enhanced) implementation of 4DD, which attempts to use the algorithm to its full potential, that is by making use of the previous velocity volume, here's some bullet points I have about it:
The following shows corrected velocities from 1040 UTC on May 20th, which I showed in a previous post. This time though I had |
Adding the capability to use a previous velocity volume as an initial condition to the dealiasing routine. Note that this enhancement was how the original 4DD algorithm was intended to be implemented, as it should provide better results. This commit attempts to address some of the issues raised in issue ARM-DOE#119.
The SW XSAPR consistently has a missing sector to its NE after |
Wonder if it is connected to the area of clutter.. I am working on some preprocessing that may help.. |
Kirk.. I noticed that.. I could play with stuff to slightly reduce it.. but always had that wedge.. I wonder if I nuke some clutter in that wedge if it would fix it… This is an interesting, if irritating, issue.. Scott Collis On Apr 15, 2014, at 2:19 PM, Kirk North notifications@github.com wrote:
|
Using only the interpolated sounding data I'm not getting a missing sector from the NW XSAPR on 2011-05-11. @kirknorth Are you using a previous volume, and did it have a missing sector? |
@jjhelmus, I was using both a previously dealiased volume (which has the same missing sector) and the interpolated sounding. And I still see remnants of that missing sector in your plot too 😉 |
I think I'll start with trying to explain this plot. The lowest sweep from the NW XSAPR (I6) using only sounding data. Reflectivity and Bergen and Albers filter turned off, yet the algorithm still masks the noon to ~1 o'clock sector. Code for my own records and if anyone is playing along at home.. import matplotlib
matplotlib.rcParams['axes.titlesize'] = 'small'
import matplotlib.pyplot as plt
import netCDF4
import pyart
SOND_NAME = 'sgpinterpolatedsondeC1.c1.20110511.000000.cdf'
RADAR_NAME = 'data/XNW110511181804.RAW7CJL'
# read in the data
radar = pyart.io.read_sigmet(RADAR_NAME)
radar = radar.extract_sweeps([0])
# find and extract sonde data
target = netCDF4.num2date(radar.time['data'][0], radar.time['units'])
interp_sounde = netCDF4.Dataset(SOND_NAME)
t = pyart.correct.find_time_in_interp_sonde(interp_sounde, target)
height, speed, direction = t
# perform dealiasing
dealias_data = pyart.correct.dealias_fourdd(
radar, sounding_heights = height * 1000.0, sounding_wind_speeds=speed,
sounding_wind_direction=direction, keep_original=False,
prep=0, filt=0)
radar.add_field('corrected_velocity', dealias_data)
# create a plots of velocity
fig = plt.figure(figsize=(10, 4))
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
display = pyart.graph.RadarDisplay(radar)
display.plot('velocity', 0, vmin=-20, vmax=20, ax=ax1,
colorbar_label='m/s', fig=fig)
display.plot('corrected_velocity', 0, vmin=-34, vmax=34, ax=ax2,
colorbar_label='m/s', mask_outside=False)
plt.show()
#fig.savefig('test.png') |
Yes definitely start with the simplest case first. During your search, can you please verify and/or check these points:
|
@kirknorth Try changing the MAXRAYS constant in FourDD.h line 20 to something like 500 and then recompile _fourdd_interface.so. From the few cases I've tried this solved the missing sector problem. It makes some sense for the XSAPRs which have 360 or 400 rays per sweep so with a limits of 350 were indexing the |
Nice work @jjhelmus. I hadn't thought to look at the MAXRAYS constant since it is not very well documented in the header file and therefore didn't really know its purpose. By the way, did you change any of the other constants to produce that latest plot? Comparing it to your previous plot (18:36 UTC) I see that not only is the NE sector filled in but also all of the other missing areas. This is promising! |
I've checked both the SW and NW X-bands and I'm quite certain that increasing MAXRAYS does get rid of our missing sector issues for both these radars. So we can tick that off the list. For the record I changed its original value of 375 to 500. That being said, it does not help with the spurious velocities issue like that showcased in the first post of this issue as shown below. I have seen these artifacts appear when using only interpolated sounding and when using both a previous volume and interpolated sounding. |
The answer to question 1 about the As the algorithms exists now Obviously this is not what should be happening during volume prep. Once I'm done with my refactor I'll be changing the code so that gates outside of a pre-defined high/low reflectivity and marked with another flag allow for marking when the reflectivity is missing. I don't want to change it yet as my regression tests that I'm using to refactor will fail. Also I'm planning of not including the |
Okay this makes some sense since when I set NODBZRMRV = 1 it does not remove gates that I've flagged as bad, i.e. with _FillValue (which then get set to the RSL bad value). Instead it only appeared to be removing gates with low reflectivity values. Now removing noisy gates is vital since the solution 4DD finds can become unstable and produce velocity values 3-4x the Nyquist within these regions and this can infiltrate cloud boundaries. Here's an example of that problem manifesting itself in gridded data: A good despeckling algorithm (like that in PR #133) can remove the isolated pixels near the radar, so I'm not worried about that. However, near the cloud boundary at (x, y) = (20, -10) in the z = 0, 1.5 km panels you can see an obvious poor correction. This is a direct result of the surrounding noisy region not being flagged as bad and its solution infiltrating the cloud edge during 4DD. This is not a result of the radius of influence during gridding capturing bad pixels within the noise, as the gates in the radial data show this as well. |
One option for filtering out bad gates until this is fixed in the algorithm is to do the masking before calling the dealias function. Something like
Would mask velocities where at "bad" reflectivity gates. You do lose the velocity data there, but if it is needed you can make a copy before masking. |
Hmm, does that even work? I removed the noise first from the reflectivity field then did is_bad = radar.fields['corrected_reflectivity']['data'].mask
radar.fields['velocity']['data'].mask = is_bad and this was the result after 4DD with sounding ICs only... You can see the original velocity field has the proper mask applied, but the corrected field in fact still retained the original velocity data which 4DD tried to correct. Leads me to believe we're not handling missing values properly. |
BTW: all who are on and following this thread.. Now is the time to send @jjhelmus your bad velocity data for testing |
@kirknorth I'm just looking at this conversation, but I've experienced some bad masking of data points. Not sure if this applies in this case, but I had a failure trying to apply the mask from one field to another as you showed in the code an hour ago. I had to use something like this to make it work: |
I wonder if this is caused by the FourDD code using an equality to test for a floating point missingValue. With roundoff and numerical precision issues using == with float is a bad idea. |
There is some literature in the numpy masked that says an equality should be tested using np.ma.masked_value (if I remember correctly). |
@jjhelmus yes one thought of mine was that there was a loss of precision between the C source code and the floating point |
Adding the capability to use a previous velocity volume as an initial condition to the dealiasing routine. Note that this enhancement was how the original 4DD algorithm was intended to be implemented, as it should provide better results. This commit attempts to address some of the issues raised in issue ARM-DOE#119.
@kirknorth You might want to have a play with my fourdd_enh branch. It is the results of refactoring the FourDD code, fixes the reflectivity filtering issue, and allow you to set a bunch of dealiasing parameters at run time to see the effects of changing them. It doesn't seem to solve the high velocities seen in the original May, 20, 2011 case, but you might be able to improve some cases. |
An IPython notebook which looks at why FourDD fails to dealias the May 20, 2011 CSAPR volume we have been looking at. Basically the algorithm is too greedy in finding gates which agree with the sounding data and if the gates agree with the sounding they are taken to be corrected even if later this is a bad choice. |
Great stuff @jjhelmus On 5/12/14, 5:24 PM, Jonathan J. Helmus wrote:
|
@kirknorth I've started to work on a Doppler dealiasing routine which uses multi-dimensional phase unwrapping for the initial (and currently only) unfolding. The technique looks to have promise, code is available in the dealias_phase_unwrap branch. Using Example of dealiasing the top most sweep from a CSAPR volume that gave FourDD issues: import pyart
import numpy as np
radar = pyart.io.read('105235.mdv')
corr_vel = pyart.correct.dealias_unwrap_phase(radar, unwrap_unit='sweep')
radar.add_field('corrected_velocity', corr_vel)
import matplotlib.pyplot as plt
display = pyart.graph.RadarDisplay(radar)
display.plot_ppi('corrected_velocity', 16, vmin=-40, vmax=40)
display.set_limits((-20, 20), (-20, 20))
plt.show() |
@jjhelmus awesome work, I'll fetch your testing branch and put it through the ringer since I've got a lot of SAPR volumes that 4DD struggled with. The May 24th case was quite a difficult case to dealias, e.g., |
Hi Kirk, Is this data from the MC3E project? I’ll be beginning some work in Sep with data from that project. Good to know the people familiar with data quirks already! On Jul 17, 2014, at 2:09 PM, Kirk North notifications@github.com wrote:
|
Oops, ignore, that was for Kirk. On Jul 17, 2014, at 2:09 PM, Kirk North notifications@github.com wrote:
|
@nguy yes sir! MC3E ran from ~ April 20th to ~ June 6th, 2011. |
@kirknorth That's great! I might be in contact with you in the future depending on which cases I use. Are you just looking at the SAPR? |
Yes, for now my focus is on the SAPR data, using it for 3D-VAR wind retrievals. |
@jjhelmus, I tested your new branch on a series of CSAPR data from May 24th, 2011 covering ~ 4 hours. I have three comments so far.
|
I agree with point 1, and I think the masking can be figured out. I didn't spend much time working out if all the masking but will in a more formalized version. I think a lot of point 2 comes from the fact the algorithm defaults to unwrapping a single sweep at a time You can override this using using Taking into account the 3D structure (or even the 4D nature if possible) should improve this greatly. There is a 3D unwrapping implementation in the branch but it seems to have some bugs with specifying which axes have periodic boundaries (causes seg faults on my machine). I addition the implementation assume a regular rectangular spacing which isn't correct for radar volumes. I have some ideas on how to address this. Never-the-less the technique looks promising. |
Two dealiasing routines have been added to Py-ART since the last update to this issue, |
Agreed.. in addition I encourage all to actively test these and let us know about On 3/3/15 9:35 AM, Jonathan J. Helmus wrote:
|
Closing this issue, feel free to reopen if further discussion is needed. |
The images below show some poor velocity unfolding behaviour. @scollis, I'm not sure if this was the same problem you were running into with the X-bands. Note how this is for the C-band.
From my quick investigation, there does appear to be some relation between areas of beam blockage and/or extinction and poor velocity correction. I'm also showing the uncorrected velocity field for reference.
In terms of the C-SAPR files that I processed for the May 20th, 2011 case, this issue was quite prevalent. In the 69 files that I processed between 0500 and 1300 UTC, about 23% had these artifacts.
The text was updated successfully, but these errors were encountered: