RadioAstronomySoftwareGroup / pyuvdata

A pythonic interface for radio astronomy interferometry data (uvfits, miriad, others)
https://pyuvdata.readthedocs.io/en/latest/index.html
BSD 2-Clause "Simplified" License
83 stars 27 forks source link

UnboundLocalError when converting a uvbeam object to healpix with nside other than default #941

Closed kbharatgehlot closed 3 years ago

kbharatgehlot commented 3 years ago
---------------------------------------------------------------------------
UnboundLocalError                         Traceback (most recent call last)
<ipython-input-7-2d0eee290741> in <module>
----> 1 uvb.to_healpix(nside=128)

~/src/anaconda/envs/HERA/lib/python3.7/site-packages/pyuvdata/uvbeam/uvbeam.py in to_healpix(self, nside, run_check, check_extra, run_check_acceptability, inplace)
   1993             run_check=run_check,
   1994             check_extra=check_extra,
-> 1995             run_check_acceptability=run_check_acceptability,
   1996         )
   1997 

~/src/anaconda/envs/HERA/lib/python3.7/site-packages/pyuvdata/uvbeam/uvbeam.py in interp(self, az_array, za_array, az_za_grid, healpix_nside, healpix_inds, freq_array, freq_interp_tol, polarizations, return_bandpass, reuse_spline, spline_opts, new_object, run_check, check_extra, run_check_acceptability)
   1804             freq_interp_kind=kind_use,
   1805             polarizations=polarizations,
-> 1806             **extra_keyword_dict,
   1807         )
   1808 

~/src/anaconda/envs/HERA/lib/python3.7/site-packages/pyuvdata/uvbeam/uvbeam.py in _interp_az_za_rect_spline(self, az_array, za_array, freq_array, freq_interp_kind, freq_interp_tol, polarizations, reuse_spline, spline_opts)
   1420 
   1421                         if do_interp:
-> 1422                             if np.iscomplexobj(data_use):
   1423                                 # interpolate real and imaginary parts separately
   1424                                 real_lut = interpolate.RectBivariateSpline(

UnboundLocalError: local variable 'data_use' referenced before assignment
bhazelton commented 3 years ago

@kbharatgehlot can you give the call that generated this error? I think I see the problem but I want to be sure and I want to add a test that breaks as a result of this bug.

kbharatgehlot commented 3 years ago

I use the following command after loading the uvbeam object: uvb.to_healpix(nside=128)

bhazelton commented 3 years ago

hmm, based on my look at the code, I suspect it has more to do with where the input beam object came from rather than the nside. Which beam object are you using?

kbharatgehlot commented 3 years ago

I see. I created this beam object from a numpy array and filled in the required parameters by myself. It did pass the check() though.

bhazelton commented 3 years ago

ah, that makes sense. Can you give me the code you used to make it? I can use that to make a test that will fail on master but should pass after my bug fix.

kbharatgehlot commented 3 years ago
beam_data = np.ones((61,181,361),dtype=float)
Nfreqs, Ntheta, Nphi = beam_data.shape

thetas = np.linspace(-0.5*np.pi,0.5*np.pi,Ntheta)
phis   = np.flip(np.linspace(0.,2*np.pi,Nphi))
freqs  = (np.arange(0.,Nfreqs) + 40.)*1.e6   

uvb = UVBeam()

uvb.Nfreqs = Nfreqs
uvb.Naxes1 = Nphi
uvb.Naxes2 = Ntheta
uvb.Nspws  = 1
uvb.Naxes_vec = 1

uvb.bandpass_array = np.ones((1,Nfreqs),dtype=float)
uvb.polarization_array = [-5]

beam_array =  np.zeros((1,1,1,Nfreqs,Ntheta,Nphi),dtype=float)
beam_array[0,0,0,:,:,:] = beam_data
uvb.data_array = beam_array
uvb.data_normalization = 'peak'
uvb.feed_name = 'Dipole'
uvb.freq_array = freqs[np.newaxis,:]
uvb.pixel_coordinate_system = 'az_za'
uvb.spw_array = [0]
uvb.Naxes1 = Nphi
uvb.Naxes2 = Ntheta
uvb.Npols = 1
uvb.axis1_array = phis
uvb.axis2_array = thetas
uvb.telescope_name = 'EDGES'

uvb.freq_interp_kind = "cubic"
uvb.interpolation_function = "az_za_simple"
uvb.x_orientation = 'north'

uvb.feed_version = '1.0'
uvb.model_name = 'FEKO-Low_newgnd'
uvb.model_version = '1.0'
uvb.history = 'FEKO simulation with new ground plane'

uvb.set_simple()
uvb.set_power()

uvb.check()

uvb.write_beamfits('test_beam.fits', run_check=True, check_extra=True, run_check_acceptability=True, clobber=True)

I replaced the beam_data array above with a generic numpy array. If I add the uvb.to_healpix(nside=128) before writing then the code crashes, however, not specifying nside parameter works normally.

kbharatgehlot commented 3 years ago

I tested the new PR, but I got the interpolation error below:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-4-2d0eee290741> in <module>
----> 1 uvb.to_healpix(nside=128)

~/src/anaconda/envs/HERA/lib/python3.7/site-packages/pyuvdata/uvbeam/uvbeam.py in to_healpix(self, nside, run_check, check_extra, run_check_acceptability, inplace)
   1997             run_check=run_check,
   1998             check_extra=check_extra,
-> 1999             run_check_acceptability=run_check_acceptability,
   2000         )
   2001 

~/src/anaconda/envs/HERA/lib/python3.7/site-packages/pyuvdata/uvbeam/uvbeam.py in interp(self, az_array, za_array, az_za_grid, healpix_nside, healpix_inds, freq_array, freq_interp_tol, polarizations, return_bandpass, reuse_spline, spline_opts, new_object, run_check, check_extra, run_check_acceptability)
   1808             freq_interp_kind=kind_use,
   1809             polarizations=polarizations,
-> 1810             **extra_keyword_dict,
   1811         )
   1812 

~/src/anaconda/envs/HERA/lib/python3.7/site-packages/pyuvdata/uvbeam/uvbeam.py in _interp_az_za_rect_spline(self, az_array, za_array, freq_array, freq_interp_kind, freq_interp_tol, polarizations, reuse_spline, spline_opts)
   1444                                     phi_use[0, :],
   1445                                     data_use[index0, index1, index2, index3, :],
-> 1446                                     **spline_opts,
   1447                                 )
   1448                                 lut = get_lambda(lut)

~/src/anaconda/envs/HERA/lib/python3.7/site-packages/scipy/interpolate/fitpack2.py in __init__(self, x, y, z, bbox, kx, ky, s)
   1167             raise ValueError('x must be strictly increasing')
   1168         if not np.all(diff(y) > 0.0):
-> 1169             raise ValueError('y must be strictly increasing')
   1170         if not ((x.min() == x[0]) and (x.max() == x[-1])):
   1171             raise ValueError('x must be strictly ascending')

ValueError: y must be strictly increasing
bhazelton commented 3 years ago

I think the problem is that you've defined the phis to run from large to small rather than small to large, which breaks scipy's interpolation. We can add a check for that kind of thing and error with a more useful error message, but I think it if you just remove the np.flip from the phi definition it will fix this error.