Sigmet/Iris KDP data is not being properly read

[johnnypeste]

In Portugal, we have several band C radars.

The data is being stored as RAW IRIS data for all radars, despite some of them are rainbow and others are vaisala.

Lisbon radar is more recent and it is originaly rainbow and then converted to iris format

Lisbon radar seems perfect when I read it in pyart with:

radar = pyart.io.read_sigmet(filename)

The kdp seems reasonable and match the results given by productx (vaisala tool)

The maximum value is for the point [30,153] with kdp = 6.71

LIS240627193002.zip

The problem is for the Porto radar which is older and is originally iris raw data.

PRT240627193006.zip

doing the same:

radar = pyart.io.read_sigmet(filename)

I get maximum kdp in the points [133,93] and [134,94] with a value of 32.706 degrees/km (!!!), which is not a reasonable value
(I think those pixels are more or less at the same place as pixel [30,153] of Lisbon radar which has 6.71 has already writen...)

Reading this Porto file with productx (vaisala tool) gives kdp = 6.14, which is similar to the 6.71 that Lisbon gives...

Now, 32.706 / 6.14 is 5.33 which the radar wavelenght (!!!)

Dividing all pyart kdp values by 5.33, I get a field which is the one decoded by Vaisala software (!!)

But only for Porto kdp data, since (as I told you) for Lisbon kdp data pyart is decoding well.

There must be some metadata information in Porto file, that pyart is desregarding, but I cannot find it...

The problem is only with kdp. Other parameters (like Z) agree in both radars and both pyart and productx

[kmuehlbauer]

Looking at the relevant code snippet, it looks like your workaround is really needed:

pyart/pyart/io/_sigmetfile.pyx

Lines 739 to 756 in 28c16de

	elif data_type_name == 'KDP':
	# KDP, 14, KDP (Differential phase) (1 byte)
	# 1-byte KDP format, section 4.3.12
	# Note that this data should be divided by the wavelength in cm
	# as is the units are deg * cm / km
	# above 128 use positive value equation
	exp = np.power(600., (ndata[ndata > 128] - 129.) / 126.)
	out[ndata > 128] = 0.25 * exp
	# below 128, use negative value equation
	exp = np.power(600., (127. - ndata[ndata < 128]) / 126.)
	out[ndata < 128] = -0.25 * exp
	# equal to 128, zero
	out[ndata == 128] = 0
	mask[ndata == 0] = True
	mask[ndata == 255] = True

LIS radar works because that KDP raw data is 2-byte, which has different decoding. PRT radar has 1-byte KDP.

[johnnypeste]

I could do that, but I need to install from source, right?

I've installed pyart with conda...

How can I change pyx file to become effective in the "import pyart" of my program?

I've tried but it does not work... I hope the devs will be quick to change the code so that conda upgrade can work

[zssherman]

@kmuehlbauer I guess i'm trying to understand when or when not to divide. Because for 2 bytes we would not want to divide? Correct me if im wrong as I'm not sure and is the part I'm still confused on. Which I was thinking if we need to do a check for when or when not to divide? Otherwise that code change you suggested would be simple enough. And @johnnypeste you would need to local install it with pip -e . etc so it cythonizes the files.

[kmuehlbauer]

@zssherman Only in the case of 1byte KDP the division by wavelength_cm is needed.

So the name is "KDP" for 1byte and "KDP2" for 2byte so it's very easy to check for.

[zssherman]

@kmuehlbauer Ahhhh that makes sense! And since we are checking for kdp and not kdp2 in that case. I can submit the change here soon. Thank you for that clarification!

[johnnypeste]

Go for it... Thanks :-) I will wait for the conda update

[johnnypeste]

yep... i need the code working for both lisbon (2 bytes) and porto (1 byte) Yahoo Mail: Search, organise, conquer On Thu, 12 Dec 2024 at 18:13, Zach ***@***.***> wrote: The code predates my knowledge of the division by wavelength. If we need to check between the two bytes, I'm trying to think of the best way to approach this. — Reply to this email directly, view it on GitHub, or unsubscribe. You are receiving this because you authored the thread.Message ID: ***@***.***>

[zssherman]

Change is in PR #1707

[johnnypeste]

Will this change be available with:

conda update -c conda-forge arm_pyart

?

When?

Thanks

[zssherman]

@johnnypeste Apologies, I'll get that out today

[zssherman]

The github release is out, just need conda to pick up the new release which can be sometime later today.

[zssherman]

And thank you @kmuehlbauer for all your help on this!

[johnnypeste]

I read the iris programmers manual. Wavelenght is stored in 1/100 of cm.
VEL and WIDTH have a factor of 10 000 in the nyquist which is correct.
When you introduced:

data['KDP'] /= wavelength_cm it should be:
data['KDP'] /= wavelength_cm/100

[zssherman]

That would make sense, do you have a test file I can work with? @johnnypeste

[johnnypeste]

changing to xradar is a bit tedious, because I have toi change the names of all variables in all my programs... I can do that, but please change the code in pyart, because it is wrong for sure...

Regards

[mgrover1]

@johnnypeste - contributions are always welcome! This is our suggested workflow moving forward - to read with xradar, then work with Py-ART.

[johnnypeste]

PRT240627193006.zip

[zssherman]

Perfect thanks!

[zssherman]

Confirmed to be the issue, adding the conversion and will do a new release.

[zssherman]

I get the matching 6.1361995 between both xradar and pyart now.

[johnnypeste]

The problem of using this code ;

import matplotlib.pyplot as plt
import xradar as xd

import pyart

read in the data

tree = xd.io.open_iris_datatree(filepath)
radar = tree.pyart.to_radar()

is that when I want to use:

pyart.map.grid_from_radars(radar,....),

I get this error:

grid = pyart.map.grid_from_radars(radar,
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/data/rna2100-A/SoftwareGCC9/anaconda3/envs/wradlib/lib/python3.11/site-packages/pyart/map/grid_mapper.py", line 86, in grid_from_radars
if len(radars) == 0:
^^^^^^^^^^^
TypeError: object of type 'Xradar' has no len()

[zssherman]

@johnnypeste Can you wrap the radar in a tuple, so its grid = pyart.map.grid_from_radars((radars,), grid_shape, etc)

[johnnypeste]

yep... that was the trick

[johnnypeste]

I talked too soon...

tree = xd.io.open_iris_datatree(filenames[0])
radar0 = tree.pyart.to_radar()

grid = pyart.map.grid_from_radars((radar0,), ...........etc

works well...

HOWEVER:

tree = xd.io.open_iris_datatree(filenames[0])
radar0 = tree.pyart.to_radar()
tree = xd.io.open_iris_datatree(filenames[1])
radar1 = tree.pyart.to_radar()

grid = pyart.map.grid_from_radars((radar0,radar1), .... etc
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/data/rna2100-A/SoftwareGCC9/anaconda3/envs/wradlib/lib/python3.11/site-packages/pyart/map/grid_mapper.py", line 93, in grid_from_radars
grids = map_gates_to_grid(radars, grid_shape, grid_limits, **kwargs)
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
File "/data/rna2100-A/SoftwareGCC9/anaconda3/envs/wradlib/lib/python3.11/site-packages/pyart/map/gates_to_grid.py", line 186, in map_gates_to_grid
radar.gate_longitude["data"], radar.gate_latitude["data"], projparams
^^^^^^^^^^^^^^^^^^^^
AttributeError: 'Xradar' object has no attribute 'gate_longitude'. Did you mean: 'gate_altitude'?

[johnnypeste]

Just to tell you that the change seems to be good now, after the division bt 100. However there are other issues which lead me to use xradar, namely the pixel position which is not correct in pyart (there is a shift of at leat 500 meters. Since you told me xradar is the guy for the job, I'm doing the shift...

Unfortunately,

grid = pyart.map.grid_from_radars((radar1,radar2), ......... etc

is not accepted when radar1 and radar2 are coming from xradar (as stated in my last post...)

and I have one sigmet file for each elevation and cannot find a way to merge them in a single file..

Attaching 3 files in a zip file
faro.zip
to see if some kind soul can help...

Regards

[zssherman]

If each radar is by sweep, can you try joining the radar by pyart.utils.join_radar. It should be able to stack them by sweep in one radar object. Grid from radars works with multiple radars but not radar objects separated by sweep i believe

[zssherman]

This is how we've merged sweep files in the past:
https://arm-development.github.io/sail-xprecip-radar/gluing_and_inventory/concatinate_file.html

[johnnypeste]

pyart.utils.join_radar doesn't work also... It complains that radar1 and radar2 (obtained with xradar) don't have sweep_number (originally they all have sweep_number = 0 but is not written when converted from xradar to pyart, like, by the way radar.gate_longitude["data"] and radar.gate_latitude["data"]

[zssherman]

for grid_from_radars

[zssherman]

what are you setting as grid_shape and grid_limits

[johnnypeste]

projection = ccrs.LambertConformal(
central_latitude=radar1.latitude["data"][0],
central_longitude=radar1.longitude["data"][0],
)

cappi_altitude_min = 1000.0
cappi_altitude_max = 1251.0

antenna_altitude = radar1.altitude["data"][0]

Z_min = cappi_altitude_min - antenna_altitude
Z_max = cappi_altitude_max - antenna_altitude

grid_resolution_z = 250.0

z_grid_limits = (Z_min , Z_max)
y_grid_limits = (-299580.9,299580.9)
x_grid_limits = (-299580.9,299580.9)

grid_resolution = 998.

def compute_number_of_points(extent, resolution):
return int((extent[1] - extent[0])/resolution)

z_grid_points = compute_number_of_points(z_grid_limits, grid_resolution_z)
x_grid_points = compute_number_of_points(x_grid_limits, grid_resolution)
y_grid_points = compute_number_of_points(y_grid_limits, grid_resolution)

grid = pyart.map.grid_from_radars((radar,1),
gridding_algo='map_gates_to_grid',
roi_func="dist_beam",
constant_roi=None,
min_radius=1000.0, h_factor=(1.0, 1.0, 1.0), nb=1.0, bsp=1.0,
grid_shape=(z_grid_points,
y_grid_points,
x_grid_points),
grid_limits=(z_grid_limits,
y_grid_limits,
x_grid_limits),
)

[mgrover1]

@johnnypeste - perhaps it would be best to set up a call to discuss this more in detail, and start opening some issues to track this?

[kmuehlbauer]

@johnnypeste FYI: The xradar IRIS reader retrieves the bin center of the range bins. The manual doesn't state if the given Range of first bin in centimeters and Range of last bin in centimeters (4.3.59 task_range_info p.60) is start/center/end of bin. But I confirmed with Vaisala that it's bin center.

See also #1329 (esp. #1329 (comment)). The implementation is here: https://github.com/openradar/xradar/blob/830d86b1c6290f1dce0e73c60a1d3b819735f906/xradar/io/backends/iris.py#L3884-L3904.

[zssherman]

@johnnypeste Heres code, let me know if it looks ok, adjust if need be. I don't know which if xradar or pyart is the correct for gate center. xradar doesn't have full functionality for all pyart functions yet, I wouldn't have suggested it if I knew we needed to join radars.

import pyart
import xradar as xd
import cartopy.crs as ccrs

# Read radars
filename1 = "FAR241115235002.RAWBX40"
filename2 = "FAR241115235023.RAWBX41"
filename3 = "FAR241115235043.RAWBX43"
radar1 = pyart.io.read(filename1)
radar2 = pyart.io.read(filename2)
radar3 = pyart.io.read(filename3)

# offset range by 500 meters to fix issue. If it does not need to be offset, remove the addition of 500.0 meters.
# Initalize the gates to fix gate spacing
radar1.range['data'] = radar1.range['data'] + 500.0
radar1.init_gate_altitude()
radar1.init_gate_longitude_latitude()
radar1.init_gate_x_y_z()

radar2.range['data'] = radar2.range['data'] + 500.0
radar2.init_gate_altitude()
radar2.init_gate_longitude_latitude()
radar2.init_gate_x_y_z()

radar3.range['data'] = radar3.range['data'] + 500.0
radar3.init_gate_altitude()
radar3.init_gate_longitude_latitude()
radar3.init_gate_x_y_z()

# combine first two sweeps
radar_mid = pyart.util.join_radar(radar1, radar2)

# add final sweep to the first two.
radar_final = pyart.util.join_radar(radar_mid, radar3)

# Shows sweeps are in order
radar_final.fixed_angle



projection = ccrs.LambertConformal(
central_latitude=radar1.latitude["data"][0],
central_longitude=radar1.longitude["data"][0],
)

cappi_altitude_min = 1000.0
cappi_altitude_max = 1251.0

antenna_altitude = radar1.altitude["data"][0]

Z_min = cappi_altitude_min - antenna_altitude
Z_max = cappi_altitude_max - antenna_altitude

grid_resolution_z = 250.0

z_grid_limits = (Z_min , Z_max)
y_grid_limits = (-299580.9,299580.9)
x_grid_limits = (-299580.9,299580.9)

grid_resolution = 998.

def compute_number_of_points(extent, resolution):
    return int((extent[1] - extent[0])/resolution)

z_grid_points = compute_number_of_points(z_grid_limits, grid_resolution_z)
x_grid_points = compute_number_of_points(x_grid_limits, grid_resolution)
y_grid_points = compute_number_of_points(y_grid_limits, grid_resolution)

grid = pyart.map.grid_from_radars((radar_final),
                                  gridding_algo='map_gates_to_grid', roi_func="dist_beam",
                                  constant_roi=None,
                                  min_radius=1000.0, h_factor=(1.0, 1.0, 1.0), nb=1.0, bsp=1.0,
                                  grid_shape=(z_grid_points, y_grid_points, x_grid_points),
                                  grid_limits=(z_grid_limits, y_grid_limits, x_grid_limits),)

display = pyart.graph.GridMapDisplay(grid)

display.plot_grid('reflectivity')

[johnnypeste]

thanks a lot... tomorrow i will check,,, But i believe that if xradar is not going to be userd, I can put firectly (radar1,radar2,radar3) in pyart.map.grid_from_radars

[zssherman]

@johnnypeste I would not, as the gridding works as if the radar1, 2 etc are different radars like if your gridding multiple radar locations in a region. Since this is one location with multiple sweeps, I would join them into one radar as I did above.

[zssherman]

Nvm I tested it, look like it did just fine, or atleast the result was the same. When plugging in each sweep. Idk use the code and adjust it to what you need

[johnnypeste]

ok... i can also check the differences,,,

pyart.io.read(filename1)

is the same as:

pyart.io.read_sigmet(filename1) ?

[zssherman]

@johnnypeste yes, so io.read is an auto reader that checks the filetype then uses the correct reader in this case read_sigmet. read_sigmet is the direct call. Either is fine.

[mgrover1]

@johnnypeste - I opened the following issue to track which components of the xradar --> py-art bridge are not working
#1713

Please feel free to update that as you find bugs.

Also a point of clarification as well - @zssherman and I are available until tomorrow. After that, we are on a holiday shutdown and will not be available until the first week of January.

[johnnypeste]

Merry Christmas and a happy new year

[johnnypeste]

I think your code does not give exactly the same values as xradar, if I use only radar1 (the only one I can use if I want to compare the location of the pixels... (the 500 m problem...)). I continue to believe more in xradar...

[johnnypeste]

Being more detailed:

In pyart, after using:

radar1.range['data'] = radar1.range['data'] + 500.0
radar1.init_gate_altitude()
radar1.init_gate_longitude_latitude()
radar1.init_gate_x_y_z()

I get the same x,y locations as with xradar

The problem is that (for the file FAR241115125502.RAWBS2J) the maximum xradar kdp is at pixel

kdp[80,25] = 16.11

and the maximum pyart kdp is at pixel

kdp[79,25] = 16.11 (one degree less in azimuth...)

I continue to think xradar is the right one and I have to find a way to use (radar1,radar2,radar3) coming from xradar