add warning when extrapolating beam

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,

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
 

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

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

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)