diff --git a/pfb/utils/beam.py b/pfb/utils/beam.py index c46071a75..34b1f4efd 100644 --- a/pfb/utils/beam.py +++ b/pfb/utils/beam.py @@ -9,7 +9,9 @@ warnings.filterwarnings("ignore", category=NumbaDeprecationWarning) from africanus.rime.fast_beam_cubes import beam_cube_dde from africanus.rime import parallactic_angles - +import pyscilog +pyscilog.init('pfb') +log = pyscilog.get_logger('BEAM') def _interp_beam_impl(freq, nx, ny, cell_deg, btype, utime=None, ant_pos=None, phase_dir=None): @@ -75,7 +77,8 @@ def _interp_beam_impl(freq, nx, ny, cell_deg, btype, beam_extents, bfreqs, lm, parangles, point_errs, ant_scale, np.array((freq,))).squeeze() - return beam_image + import ipdb; ipdb.set_trace() + return beam_image.squeeze() def interp_beam(freq, nx, ny, cell_deg, btype, @@ -115,8 +118,15 @@ def interp_beam(freq, nx, ny, cell_deg, btype, def _eval_beam(beam_image, l_in, m_in, l_out, m_out): - beamo = RGI((l_in, m_in), beam_image, - bounds_error=True, method='linear') + try: + beamo = RGI((l_in, m_in), beam_image, + bounds_error=True, method='linear') + except: + print("Bounds error raised in beam evaluation. " + "Consider setting init.max_field_of_view > grid.field_of_view", + file=log) + beamo = RGI((l_in, m_in), beam_image, + bounds_error=False, method='linear', fill_value=None) if l_out.ndim == 2: ll = l_out mm = m_out @@ -129,7 +139,8 @@ def _eval_beam(beam_image, l_in, m_in, l_out, m_out): def eval_beam(beam_image, l_in, m_in, l_out, m_out): - nxo, nyo = l_out.shape + nxo = l_out.shape[0] + nyo = m_out.shape[-1] return da.blockwise(_eval_beam, 'xy', beam_image, 'xy', l_in, None, diff --git a/pfb/utils/misc.py b/pfb/utils/misc.py index a9e71a4b6..5e2656c7a 100644 --- a/pfb/utils/misc.py +++ b/pfb/utils/misc.py @@ -885,13 +885,6 @@ def concat_chan(xds, nband_out=1): time_max.append(ds.time_max) time_min.append(ds.time_min) - - # LB - we should be able to avoid this stack operation by using Jon's *() magic - # wgt = da.stack([ds.WEIGHT.data for ds in xdst]).rechunk(-1, -1) # verify chunking over freq axis - # vis = da.stack([ds.VIS.data for ds in xdst]).rechunk(-1, -1) - # mask = da.stack([ds.MASK.data for ds in xdst]).rechunk(-1, -1) - # freq = da.stack([ds.FREQ.data for ds in xdst]).rechunk(-1) - nrow = xdst[0].row.size nchan = freqsb.size diff --git a/pfb/workers/grid.py b/pfb/workers/grid.py index 2e49875dc..cbea26e85 100644 --- a/pfb/workers/grid.py +++ b/pfb/workers/grid.py @@ -382,10 +382,10 @@ def _grid(xdsi=None, **kw): cell_deg = np.rad2deg(cell_rad) l = (-(nx//2) + da.arange(nx)) * cell_deg + np.deg2rad(x0) m = (-(ny//2) + da.arange(ny)) * cell_deg + np.deg2rad(y0) - ll, mm = da.meshgrid(l, m, indexing='ij') + # ll, mm = da.meshgrid(l, m, indexing='ij') l_beam = ds.l_beam.data m_beam = ds.m_beam.data - bvals = eval_beam(ds.BEAM.data, l_beam, m_beam, ll, mm) + bvals = eval_beam(ds.BEAM.data, l_beam, m_beam, l, m) out_ds = out_ds.assign(**{'BEAM': (('x', 'y'), bvals)}) # get the model diff --git a/tests/test_clean.py b/tests/test_clean.py index ed76cf4e8..e274e08cd 100644 --- a/tests/test_clean.py +++ b/tests/test_clean.py @@ -189,7 +189,7 @@ def test_clean(do_gains, ms_name): init_args["data_column"] = "DATA" init_args["flag_column"] = 'FLAG' init_args["gain_table"] = gain_path - init_args["max_field_of_view"] = fov + init_args["max_field_of_view"] = fov*1.1 init_args["overwrite"] = True init_args["channels_per_image"] = 1 from pfb.workers.init import _init diff --git a/tests/test_spotless.py b/tests/test_spotless.py index 43bf634ef..12d4a1199 100644 --- a/tests/test_spotless.py +++ b/tests/test_spotless.py @@ -293,7 +293,7 @@ def test_spotless(ms_name): # init_args["weight_column"] = 'WEIGHT_SPECTRUM' init_args["flag_column"] = 'FLAG' init_args["gain_table"] = None - init_args["max_field_of_view"] = fov + init_args["max_field_of_view"] = fov*1.1 init_args["overwrite"] = True init_args["channels_per_image"] = 1 xdso2 = _init(**init_args)