Python 3.6 bug in gridding routines

[tjlang]

Gridding routines appear to be broken in Python 3.6. In order to account for mobile radars, Py-ART Radar object allows arrays for radar lat/lon/alt, altho a lot of code heritage assumes scalars for these parameters. It appears that doing a float conversion on the z_disp calculation, or if that is removed, the appending immediately after this calculation raises a TypeError in 3.6 due to this array/scalar confusion. Here is the error below. This completely breaks gridding.

---

TypeError Traceback (most recent call last)

in () 1 g1 = grid_radar(r1, origin=(r2.latitude['data'][0], r2.longitude['data'][0]), ----> 2 xlim=(-100000, 50000), ylim=(0, 150000), grid_shape=(20, 151, 151)) 3 g2 = grid_radar(r2, origin=(r2.latitude['data'][0], r2.longitude['data'][0]), 4 xlim=(-100000, 50000), ylim=(0, 150000), grid_shape=(20, 151, 151)) in grid_radar(radar, grid_shape, xlim, ylim, zlim, fields, origin) 11 grid_limits=(zlim, ylim, xlim), 12 grid_origin=origin, fields=fields, ---> 13 gridding_algo='map_gates_to_grid', grid_origin_alt=0.0) 14 print(time.time()-bt, 'seconds to grid radar') 15 return grid /Users/tjlang/anaconda/envs/py3/lib/python3.6/site-packages/pyart/map/grid_mapper.py in grid_from_radars(radars, grid_shape, grid_limits, gridding_algo, **kwargs) 85 grids = map_to_grid(radars, grid_shape, grid_limits, **kwargs) 86 elif gridding_algo == 'map_gates_to_grid': ---> 87 grids = map_gates_to_grid(radars, grid_shape, grid_limits, **kwargs) 88 else: 89 raise ValueError('invalid gridding_algo') /Users/tjlang/anaconda/envs/py3/lib/python3.6/site-packages/pyart/map/gates_to_grid.py in map_gates_to_grid(radars, grid_shape, grid_limits, grid_origin, grid_origin_alt, grid_projection, fields, gatefilters, map_roi, weighting_function, toa, roi_func, constant_roi, z_factor, xy_factor, min_radius, h_factor, nb, bsp, **kwargs) 97 fields = _determine_fields(fields, radars) 98 grid_starts, grid_steps = _find_grid_params(grid_shape, grid_limits) ---> 99 offsets = _find_offsets(radars, projparams, grid_origin_alt) 100 roi_func = _parse_roi_func(roi_func, constant_roi, z_factor, xy_factor, 101 min_radius, h_factor, nb, bsp, offsets) /Users/tjlang/anaconda/envs/py3/lib/python3.6/site-packages/pyart/map/gates_to_grid.py in _find_offsets(radars, projparams, grid_origin_alt) 223 z_disp = np.array(radar.altitude['data'], dtype='float') - \ 224 grid_origin_alt --> 225 offsets.append((z_disp, float(y_disp), float(x_disp))) 226 return offsets 227 TypeError: only length-1 arrays can be converted to Python scalars

[tjlang]

Looking at this in more detail, it looks like the issue will ultimately come down to fixing the Cython code. If I successfully convert the x/y/z_disp arrays coming out of _find_offsets (by just multiplying by 1.0 to ensure they are floats), I still get hamstrung by the following error.

Also, this is not just a 3.6 issue. This is also cropping up in 2.7. It appears to be related to whether or not radar lat/lon/alt each have one or many elements. If many, gridding fails.

---

TypeError Traceback (most recent call last)

in () 1 g1 = grid_radar(r1, origin=(r2.latitude['data'][0], r2.longitude['data'][0]), ----> 2 xlim=(-100000, 50000), ylim=(0, 150000), grid_shape=(20, 151, 151)) 3 g2 = grid_radar(r2, origin=(r2.latitude['data'][0], r2.longitude['data'][0]), 4 xlim=(-100000, 50000), ylim=(0, 150000), grid_shape=(20, 151, 151)) in grid_radar(radar, grid_shape, xlim, ylim, zlim, fields, origin) 11 grid_limits=(zlim, ylim, xlim), 12 grid_origin=origin, fields=fields, ---> 13 gridding_algo='map_gates_to_grid', grid_origin_alt=0.0) 14 print(time.time()-bt, 'seconds to grid radar') 15 return grid /Users/tjlang/anaconda/lib/python2.7/site-packages/pyart/map/grid_mapper.pyc in grid_from_radars(radars, grid_shape, grid_limits, gridding_algo, **kwargs) 85 grids = map_to_grid(radars, grid_shape, grid_limits, **kwargs) 86 elif gridding_algo == 'map_gates_to_grid': ---> 87 grids = map_gates_to_grid(radars, grid_shape, grid_limits, **kwargs) 88 else: 89 raise ValueError('invalid gridding_algo') /Users/tjlang/anaconda/lib/python2.7/site-packages/pyart/map/gates_to_grid.pyc in map_gates_to_grid(radars, grid_shape, grid_limits, grid_origin, grid_origin_alt, grid_projection, fields, gatefilters, map_roi, weighting_function, toa, roi_func, constant_roi, z_factor, xy_factor, min_radius, h_factor, nb, bsp, **kwargs) 99 offsets = _find_offsets(radars, projparams, grid_origin_alt) 100 roi_func = _parse_roi_func(roi_func, constant_roi, z_factor, xy_factor, --> 101 min_radius, h_factor, nb, bsp, offsets) 102 103 # prepare grid storage arrays /Users/tjlang/anaconda/lib/python2.7/site-packages/pyart/map/gates_to_grid.pyc in _parse_roi_func(roi_func, constant_roi, z_factor, xy_factor, min_radius, h_factor, nb, bsp, offsets) 261 roi_func = DistRoI(z_factor, xy_factor, min_radius, offsets) 262 elif roi_func == 'dist_beam': --> 263 roi_func = DistBeamRoI(h_factor, nb, bsp, min_radius, offsets) 264 else: 265 raise ValueError('unknown roi_func: %s' % roi_func) pyart/map/_gate_to_grid_map.pyx in pyart.map._gate_to_grid_map.DistBeamRoI.__init__ (pyart/map/_gate_to_grid_map.c:3229)() TypeError: only length-1 arrays can be converted to Python scalars

[tjlang]

Pinging Py-ART folks - have you been able to reproduce this issue, and/or have a plan to address it?

[zssherman]

@tjlang I will try to reproduce the issue, I'll also talk to @scollis and discuss a way to address it. Thanks for bringing this up and I apologize for the delayed response.

[scollis]

Thanks Tim.. Since I have not had issues I missed this

Scott Collis

On Apr 18, 2018 at 11:42 AM, Zach Sherman notifications@github.com wrote:

@tjlang https://github.com/tjlang I will try to reproduce the issue, I'll also talk to @scollis https://github.com/scollis and discuss a way to address it. Thanks for bringing this up and I apologize for the delayed response.

— You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub https://github.com/ARM-DOE/pyart/issues/685#issuecomment-382451970, or mute the thread https://github.com/notifications/unsubscribe-auth/AAyYB4AReVDBZ1m0fbKK1z_I8j4Ce-i3ks5tp20RgaJpZM4PXkXs .

[zssherman]

@tjlang I tried to reproduce the error in 2.7 and 3.6 but did not get an error. I used the files in the multidop notebook before the radx changes to nc, 20090410_173943_ARMOR_v001_PPI.uf.gz and 20090410_174143_KHTX_v537_SUR.uf.gz. I tried changing the float conversion of z_disp or removing it, but the grid code still ran. I wonder if it could possibly be a dependency issue.

[dcon691]

I thought the problem might have been arising from using Radx. After editing sweeps using solo3, I tried converting the sweeps into a volume using Radx to a UF and cfradial file. Both of those yield the TypeError: only length-1 arrays can be converted to Python scalars. Using the raw data from NCDC, it works fine and I don't get that issue.

However, this isn't an issue when following the same procedure for the ARMOR radar. That reads and grids just fine using the UF file for the Barnes scheme.

Here is a link to the converted UF files that I've been using. https://drive.google.com/file/d/1gQLAA70ghoy2VOcA2FomzvEn8te3X-VG/view?usp=sharing

[scollis]

Ok.. Let me wrap my head around this.. @tjlang the issue is with files that were gridded with RadX?

[tjlang]

The issue I had, was that reading the UF files after conversion from RadX worked fine. Reading the UF files straight caused the Py-ART gridding issues. So kind of the opposite issue that @dcon691 had.

[scollis]

ok.. So it may be an issue in the UF reader. Can you put the file somewhere?

[tjlang]

UF files and CF radial files obtained via RadX can be found here: https://github.com/nasa/MultiDop/tree/master/examples

[scollis]

Ok.. I am not getting an error with 20090410_173943_ARMOR_v001_PPI.uf

import pyart
from matplotlib import pyplot as plt
import numpy as np
import time
%matplotlib inline

r2 = pyart.io.read('/Users/scollis/Downloads/20090410_173943_ARMOR_v001_PPI.uf')

def grid_radar(radar, grid_shape=(20, 301, 301), xlim=(-150000, 150000),
               ylim=(-150000, 150000), zlim=(1000, 20000),
               fields=['reflectivity'], origin=None):
    bt = time.time()
    radar_list = [radar]
    if origin is None:
        origin = (radar.latitude['data'][0],
                  radar.longitude['data'][0])
    grid = pyart.map.grid_from_radars(
        radar_list, grid_shape=grid_shape,
        grid_limits=(zlim, ylim, xlim),
        grid_origin=origin, fields=fields,
        gridding_algo='map_gates_to_grid', grid_origin_alt=0.0)
    print(time.time()-bt, 'seconds to grid radar')
    return grid

g2 = grid_radar(r2, origin=(r2.latitude['data'][0], r2.longitude['data'][0]),
                xlim=(-100000, 50000), ylim=(0, 150000), grid_shape=(20, 151, 151))

Returns 2.334674835205078 seconds to grid radar

My environment used is here: https://github.com/EVS-ATMOS/environments/blob/master/radar/pyart_full.yml

[tjlang]

LOL, now I can't even read the UF directly using pyart.io.read. So I can't even try to reproduce my old gridding error.

Honestly, I have never liked UF ingest in Py-ART. It just is not bulletproof. It works sometimes and not others, and how well it works changes over time as packages are gradually updated, and it even changes from UF volume to UF volume. The read_radx option makes such a difference there, and for me is close to 100% reliable. I can grid no problem with UF data ingested that way.

In terms of this particular thread's issue, I dunno what to do. Maybe leave it open and see if others come to it with similar problems, and close it in a few months if not? Ultimately, I think the issue is related to UF, that Py-ART is not very stable in terms of its ability to handle UF (i.e., it appears to be environment and architecture specific - people can't reliably reproduce each other's errors), and how much time do you want to invest in getting Py-ART to read and grid decades-old legacy radar formats for all users, when read_radx is the obvious alternative?

Ingest error below, for giggles:

---

ValueError Traceback (most recent call last)

in () ----> 1 r1 = pyart.io.read(f1[0]) 2 r2 = pyart.io.read(f2[0]) 3 print(r1.fields.keys(), r2.fields.keys()) ~/anaconda/envs/py3/lib/python3.6/site-packages/pyart/io/auto_read.py in read(filename, use_rsl, **kwargs) 137 return read_rsl(filename, **kwargs) 138 else: --> 139 return read_uf(filename, **kwargs) 140 141 # RSL only supported file formats ~/anaconda/envs/py3/lib/python3.6/site-packages/pyart/io/uf.py in read_uf(filename, field_names, additional_metadata, file_field_names, exclude_fields, delay_field_loading, **kwargs) 173 continue 174 field_dic = filemetadata(field_name) --> 175 field_dic['data'] = ufile.get_field_data(uf_field_number) 176 field_dic['_FillValue'] = get_fillvalue() 177 fields[field_name] = field_dic ~/anaconda/envs/py3/lib/python3.6/site-packages/pyart/io/uffile.py in get_field_data(self, field_number) 197 ray_data = ray.field_raw_data[field_number] 198 bins = len(ray_data) --> 199 raw_data[i, :bins] = ray.field_raw_data[field_number] 200 raw_data[i, bins:] = missing_data_value 201 ValueError: could not broadcast input array from shape (1832) into shape (1192)

[scollis]

Was this the same file I was testing? Looks like it has a variable number of gates per field or suchlike.

[tjlang]

It was the KHTX file that failed on ingest, not the ARMOR one. Obviously, we know we can read NEXRAD that hasn't been converted to UF first just fine. So now I wonder if this isn't an artifact of that. I got the file from someone else who had done processing on the volume and saved it as a UF. NEXRAD has funky enough scanning and file structure that I can see UF struggling with that. This might explain the issues @dcon691 also had with the NEXRAD file, yet everyone (including me) seems to not have problems with the ARMOR file.

Py-ART's interface to UF has always been a little touchy, so if you take a UF-converted NEXRAD file then outright failure or at least instability and platform-dependence on ingest (or gridding) makes sense. The obvious solution is first convert the file or go grab the original KHTX from AWS/NCEI. What I did for the MultiDop example was also provide CF radial versions of the volumes.

So the issue may be more NEXRAD/UF incompatibility rather than Py-ART per se. You are welcome to close this issue unless you really want to delve into the details of NEXRAD to UF conversion, which I don't think anyone has worried about in a few years.

[mjskier]

Does pyart support gridding of airborn data?

I have a cfradial file that produces the same TypeError: Only length-1 arrays can be converted to Python scalars when I call pyart.map.grid_from_radars.

If I don't set grid_origin_alt and grid_origin, it fails trying to convert a numpy array in map_gates_to_grid at line 92:

    if grid_origin_alt is None:
        grid_origin_alt = float(radars[0].altitude['data'])

If I set grid_origin_alt and grid_origin (as in Scott's May 22 answer above) if fails later in _find_offsets (gates_to_grid.py line 223)

       z_disp = float(radar.altitude['data']) - grid_origin_alt

/usr/local/anaconda3/bin/conda list | grep pyart
arm_pyart                 1.9.2                    py35_0    conda-forge

[scollis]

@nguy can you help?