@@ -321,7 +321,7 @@ def __str__(self):
         s += "\tDirection: %s\n" % self.direction
         s += "\tShape: %s\n" % self.shape
         s += "\tParams: %s\n" % self.params
-        s += "\tPolarisation frame: %s\n" % self.polarisation_frame.type
+        s += "\tPolarisation frame: %s\n" % str(self.polarisation_frame.type)
         return s
 
 
 
@@ -125,8 +125,8 @@ def deconvolve_cube(dirty: Image, psf: Image, **kwargs) -> (Image, Image):
                 else:
                     log.info("deconvolve_cube: Skipping pol %d, channel %d" % (pol, channel))
 
-        comp_image = create_image_from_array(comp_array, dirty.wcs)
-        residual_image = create_image_from_array(residual_array, dirty.wcs)
+        comp_image = create_image_from_array(comp_array, dirty.wcs, dirty.polarisation_frame)
+        residual_image = create_image_from_array(residual_array, dirty.wcs, dirty.polarisation_frame)
 
     elif algorithm == 'msmfsclean' or algorithm == 'mfsmsclean' or algorithm == 'mmclean':
         findpeak = get_parameter(kwargs, "findpeak", 'ARL')
@@ -172,8 +172,8 @@ def deconvolve_cube(dirty: Image, psf: Image, **kwargs) -> (Image, Image):
             else:
                 log.info("deconvolve_cube: Skipping pol %d" % (pol))
 
-        comp_image = create_image_from_array(comp_array, dirty_taylor.wcs)
-        residual_image = create_image_from_array(residual_array, dirty_taylor.wcs)
+        comp_image = create_image_from_array(comp_array, dirty_taylor.wcs, dirty.polarisation_frame)
+        residual_image = create_image_from_array(residual_array, dirty_taylor.wcs, dirty.polarisation_frame)
 
         return_moments = get_parameter(kwargs, "return_moments", False)
         if not return_moments:
@@ -211,8 +211,8 @@ def deconvolve_cube(dirty: Image, psf: Image, **kwargs) -> (Image, Image):
                 else:
                     log.info("deconvolve_cube: Skipping pol %d, channel %d" % (pol, channel))
 
-        comp_image = create_image_from_array(comp_array, dirty.wcs)
-        residual_image = create_image_from_array(residual_array, dirty.wcs)
+        comp_image = create_image_from_array(comp_array, dirty.wcs, dirty.polarisation_frame)
+        residual_image = create_image_from_array(residual_array, dirty.wcs, dirty.polarisation_frame)
     else:
         raise ValueError('deconvolve_cube: Unknown algorithm %s' % algorithm)
 
 
@@ -63,7 +63,7 @@ def image_raster_iter(im: Image, **kwargs) -> Image:
             wcs.wcs.crpix[1] -= y
 
             # Yield image from slice (reference!)
-            yield create_image_from_array(im.data[..., y:y + dy, x:x + dx], wcs)
+            yield create_image_from_array(im.data[..., y:y + dy, x:x + dx], wcs, im.polarisation_frame)
 
 
 def image_channel_iter(im: Image, subimages=1) -> Image:
@@ -100,4 +100,4 @@ def image_channel_iter(im: Image, subimages=1) -> Image:
         wcs.wcs.crpix[3] -= channel
 
         # Yield image from slice (reference!)
-        yield create_image_from_array(im.data[channel:channel_max, ...], wcs)
+        yield create_image_from_array(im.data[channel:channel_max, ...], wcs, im.polarisation_frame)
@@ -33,8 +33,7 @@ def image_sizeof(im: Image):
     return im.size()
 
 
-def create_image_from_array(data: numpy.array, wcs: WCS = None,
-                            polarisation_frame=PolarisationFrame('stokesI')) -> Image:
+def create_image_from_array(data: numpy.array, wcs: WCS, polarisation_frame: PolarisationFrame) -> Image:
     """ Create an image from an array and optional wcs
 
     :param data: Numpy.array
@@ -211,7 +210,8 @@ def reproject_image(im: Image, newwcs: WCS, shape=None) -> (Image, Image):
     assert isinstance(im, Image), im
     rep, foot = reproject_interp((im.data, im.wcs), newwcs, shape, order='bicubic',
                                  independent_celestial_slices=True)
-    return create_image_from_array(rep, newwcs), create_image_from_array(foot, newwcs)
+    return create_image_from_array(rep, newwcs, im.polarisation_frame), create_image_from_array(foot, newwcs,
+                                                                                                im.polarisation_frame)
 
 
 def checkwcs(wcs1, wcs2):
@@ -513,7 +513,7 @@ def fft_image(im, template_image=None):
         ft_wcs.wcs.cdelt[0] = 1.0 / (ft_shape[3] * d2r * im.wcs.wcs.cdelt[0])
         ft_wcs.wcs.cdelt[1] = 1.0 / (ft_shape[2] * d2r * im.wcs.wcs.cdelt[1])
         ft_data = ifft(im.data.astype('complex'))
-        return create_image_from_array(ft_data, wcs=ft_wcs)
+        return create_image_from_array(ft_data, wcs=ft_wcs, polarisation_frame=im.polarisation_frame)
     elif im.wcs.wcs.ctype[0] == 'UU' and im.wcs.wcs.ctype[1] == 'VV':
         ft_wcs.wcs.crval[0] = template_image.wcs.wcs.crval[0]
         ft_wcs.wcs.crval[1] = template_image.wcs.wcs.crval[1]
@@ -524,7 +524,7 @@ def fft_image(im, template_image=None):
         ft_wcs.wcs.cdelt[0] = template_image.wcs.wcs.cdelt[0]
         ft_wcs.wcs.cdelt[1] = template_image.wcs.wcs.cdelt[1]
         ft_data = fft(im.data.astype('complex'))
-        return create_image_from_array(ft_data, wcs=ft_wcs)
+        return create_image_from_array(ft_data, wcs=ft_wcs, polarisation_frame=im.polarisation_frame)
     else:
         raise NotImplementedError("Cannot FFT specified axes")
 
@@ -622,7 +622,7 @@ def convert_image_to_kernel(im: Image, oversampling, kernelwidth):
                     slicex = slice(xend + x, xstart + x, -oversampling)
                     newdata[chan, pol, y, x, ...] = im.data[chan, pol, slicey, slicex]
 
-    return create_image_from_array(newdata, newwcs)
+    return create_image_from_array(newdata, newwcs, polarisation_frame=im.polarisation_frame)
 
 
 def create_w_term_like(im: Image, w, phasecentre=None, remove_shift=False, dopol=False) -> Image:
 
@@ -63,6 +63,7 @@ def shift_vis_to_image(vis: Visibility, im: Image, tangent: bool = True, inverse
     # frames is defined.
 
     image_phasecentre = pixel_to_skycoord(nx // 2 + 1, ny // 2 + 1, im.wcs, origin=1)
+    print("Image phase centre is %s" % str(image_phasecentre))
 
     if vis.phasecentre.separation(image_phasecentre).rad > 1e-15:
         if inverse:
@@ -205,15 +206,15 @@ def invert_2d_base(vis: Visibility, im: Image, dopsf: bool = False, normalize: b
     if imaginary:
         log.debug("invert_2d_base: retaining imaginary part of dirty image")
         result = extract_mid(ifft(imgridpad) * gcf, npixel=nx)
-        resultreal = create_image_from_array(result.real, im.wcs)
-        resultimag = create_image_from_array(result.imag, im.wcs)
+        resultreal = create_image_from_array(result.real, im.wcs, im.polarisation_frame)
+        resultimag = create_image_from_array(result.imag, im.wcs, im.polarisation_frame)
         if normalize:
             resultreal = normalize_sumwt(resultreal, sumwt)
             resultimag = normalize_sumwt(resultimag, sumwt)
         return resultreal, sumwt, resultimag
     else:
         result = extract_mid(numpy.real(ifft(imgridpad)) * gcf, npixel=nx)
-        resultimage = create_image_from_array(result, im.wcs)
+        resultimage = create_image_from_array(result, im.wcs, im.polarisation_frame)
         if normalize:
             resultimage = normalize_sumwt(resultimage, sumwt)
         return resultimage, sumwt
 
@@ -341,10 +341,10 @@ def pwr23(n):
     log.info('advice_wide_field: Number of planes in w stack %d' % (vis_slices))
     log.info('advice_wide_field: Number of planes in w projection %d' % wprojection_planes)
     if wprojection_planes > 1:
-        log.info('advice_wide_field: Recommend wprojection gridding')
+        log.info('advice_wide_field: Recommend that wprojection gridding is used')
         kernel = 'wprojection'
     else:
-        log.info('advice_wide_field: Recommend 2d gridding (i.e. no wprojection)')
+        log.info('advice_wide_field: Recommend that 2d gridding (i.e. no wprojection) is used')
         kernel = '2d'
 
     return locals()
 
@@ -147,7 +147,7 @@ def comp_prop(comp, prop_name):
     return comps
 
 
-def apply_beam_to_skycomponent(sc: Union[Skycomponent, List[Skycomponent]], beam: Image, flux_limit=1.0) \
+def apply_beam_to_skycomponent(sc: Union[Skycomponent, List[Skycomponent]], beam: Image, flux_limit=0.0) \
         -> Union[Skycomponent, List[Skycomponent]]:
     """ Insert a Skycomponet into an image
     :param beam:
 
@@ -50,6 +50,7 @@
 from arl.util.coordinate_support import xyz_at_latitude
 from arl.visibility.base import create_blockvisibility
 from arl.visibility.coalesce import convert_visibility_to_blockvisibility
+from arl.skycomponent.operations import insert_skycomponent
 
 log = logging.getLogger(__name__)
 
@@ -349,7 +350,8 @@ def create_low_test_image_composite(npixel=16384, polarisation_frame=Polarisatio
 
 def create_low_test_image_from_gleam(npixel=512, polarisation_frame=PolarisationFrame("stokesI"), cellsize=0.000015,
                                      frequency=numpy.array([1e8]), channel_bandwidth=numpy.array([1e6]),
-                                     phasecentre=None, kind='cubic', applybeam=False, flux_limit=0.1) -> Image:
+                                     phasecentre=None, kind='cubic', applybeam=False, flux_limit=0.1,
+                                     radius=None, insert_method='Nearest') -> Image:
     """Create LOW test image from the GLEAM survey
 
     Stokes I is estimated from a cubic spline fit to the measured fluxes. The polarised flux is always zero.
@@ -371,25 +373,21 @@ def create_low_test_image_from_gleam(npixel=512, polarisation_frame=Polarisation
     :return: Image
     
     """
-    
-    
-    fitsfile = arl_path("data/models/GLEAM_EGC.fits")
-    
-    hdulist = fits.open(fitsfile, lazy_load_hdus=False)
-    recs = hdulist[1].data[0].array
-    ras = recs['RAJ2000']
-    decs = recs['DEJ2000']
-    
+
     if phasecentre is None:
         phasecentre = SkyCoord(ra=+15.0 * u.deg, dec=-35.0 * u.deg, frame='icrs', equinox='J2000')
-    
+        
+    if radius is None:
+        radius = npixel * cellsize / numpy.sqrt(2.0)
+
+    sc = create_low_test_skycomponents_from_gleam(flux_limit=flux_limit, polarisation_frame=polarisation_frame,
+                                                  frequency=frequency, phasecentre=phasecentre,
+                                                  kind=kind, radius=radius)
     if polarisation_frame is None:
         polarisation_frame = PolarisationFrame("stokesI")
 
     npol = polarisation_frame.npol
-    
     nchan = len(frequency)
-    
     shape = [nchan, npol, npixel, npixel]
     w = WCS(naxis=4)
     # The negation in the longitude is needed by definition of RA, DEC
@@ -398,59 +396,18 @@ def create_low_test_image_from_gleam(npixel=512, polarisation_frame=Polarisation
     w.wcs.ctype = ["RA---SIN", "DEC--SIN", 'STOKES', 'FREQ']
     w.wcs.crval = [phasecentre.ra.deg, phasecentre.dec.deg, 1.0, frequency[0]]
     w.naxis = 4
-    
     w.wcs.radesys = 'ICRS'
     w.wcs.equinox = 2000.0
 
     model = create_image_from_array(numpy.zeros(shape), w, polarisation_frame=polarisation_frame)
 
-    p = w.sub(2).wcs_world2pix(numpy.array(ras), numpy.array(decs), 1)
-    
-    p = numpy.array(p)
-    mask = numpy.isfinite(p[0, :])
-    p = [p[0, mask], p[1, mask]]
-    
-    ip = numpy.round(p).astype('int')
-    ok = numpy.where((0 <= ip[0, :]) & (npixel > ip[0, :]) & (0 <= ip[1, :]) & (npixel > ip[1, :]))[0]
-    ps = ip[:, ok]
-    
-    # For every source, we read all measured fluxes and interpolate to the
-    # required frequencies
-    fluxes = []
-    gleam_freqs = numpy.array([76, 84, 92, 99, 107, 115, 122, 130, 143, 151, 158, 166, 174, 181, 189, 197, 204,
-                               212, 220, 227])
-    gleam_flux_freq = numpy.zeros([len(ok), len(gleam_freqs)])
-    for i, f in enumerate(gleam_freqs):
-        gleam_flux_freq[:,i] = recs['int_flux_%03d' % (f)][:]
-
+    model = insert_skycomponent(model, sc, insert_method=insert_method)
     if applybeam:
         beam = create_low_test_beam(model)
+        model.data[...] *= beam.data[...]
 
-    for source in ok:
+    log.info(qa_image(model, context='create_low_test_image_from_gleam'))
 
-        fint = interpolate.interp1d(gleam_freqs * 1.0e6, gleam_flux_freq[source,:], kind=kind)
-        factual = fint(frequency)
-        if applybeam:
-            factual *= beam.data[:, 0, ip[1, source], ip[0, source]]
-        
-        fluxes.append(factual)
-    
-    fluxes = numpy.array(fluxes)
-    actual_flux = numpy.sum(fluxes)
-    
-    log.info('create_low_test_image_from_gleam: %d sources inside the image' % (ps.shape[1]))
-    
-    log.info(
-        'create_low_test_image_from_gleam: Average flux per channel in image = %.3f' % (actual_flux / float(nchan)))
-    for iflux, flux in enumerate(fluxes):
-        if not numpy.isnan(flux).any() and flux.any() > flux_limit:
-            model.data[:, 0, ps[1, iflux], ps[0, iflux]] = flux[:]
-                
-    log.info(qa_image(model), context='create_low_test_image_from_gleam')
-    print(qa_image(model), context='create_low_test_image_from_gleam')
-
-    hdulist.close()
-    
     return model
 
 
@@ -470,6 +427,7 @@ def create_low_test_skycomponents_from_gleam(flux_limit=0.1, polarisation_frame=
     catalogue. Hurley-Walker N., et al., Mon. Not. R. Astron. Soc., 464, 1146-1167 (2017), 2017MNRAS.464.1146H
 
 
+    :rtype: Union[None, List[arl.data.data_models.Skycomponent], List]
     :param flux_limit: Only write components brighter than this (Jy)
     :param polarisation_frame: Polarisation frame (default PolarisationFrame("stokesI"))
     :param frequency: Frequencies at which the flux will be estimated
@@ -481,27 +439,26 @@ def create_low_test_skycomponents_from_gleam(flux_limit=0.1, polarisation_frame=
 
     fitsfile = arl_path("data/models/GLEAM_EGC.fits")
 
-    decmin = phasecentre.dec.to('rad').value - radius/2.0
-    decmax = phasecentre.dec.to('rad').value + radius/2.0
+    rad2deg = 180.0/numpy.pi
+    decmin = phasecentre.dec.to('deg').value - rad2deg * radius/2.0
+    decmax = phasecentre.dec.to('deg').value + rad2deg * radius/2.0
 
     hdulist = fits.open(fitsfile, lazy_load_hdus=False)
     recs = hdulist[1].data[0].array
-    print("Before filter %s" % recs.shape)
+
+    # Do the simple forms of filtering in pyfits. Filtering on radious is done below.
     fluxes = recs['peak_flux_wide']
 
     mask = fluxes > flux_limit
     filtered_recs = recs[mask]
-    print("After flux filter %s" % filtered_recs.shape)
 
     decs = filtered_recs['DEJ2000']
     mask = decs > decmin
     filtered_recs = filtered_recs[mask]
-    print("After dec min filter %s" % filtered_recs.shape)
 
     decs = filtered_recs['DEJ2000']
     mask = decs < decmax
     filtered_recs = filtered_recs[mask]
-    print("After dec max filter %s" % filtered_recs.shape)
 
     ras = filtered_recs['RAJ2000']
     decs = filtered_recs['DEJ2000']
@@ -579,7 +536,7 @@ def create_low_test_beam(model: Image) -> Image:
         beam2dwcs.wcs.ctype = model.wcs.sub(2).wcs.ctype
         model2dwcs.wcs.crpix = [model.shape[2] // 2 + 1, model.shape[3] // 2 + 1]
 
-        beam2d = create_image_from_array(beam.data[0, 0, :, :], beam2dwcs)
+        beam2d = create_image_from_array(beam.data[0, 0, :, :], beam2dwcs, model.polarisation_frame)
         reprojected_beam2d, footprint = reproject_image(beam2d, model2dwcs, shape=model2dshape)
         assert numpy.max(footprint.data) > 0.0, "No overlap between beam and model"