Data-simulation, difference-imaging and sourcefinding

A tweakable end-to-end runthrough



In [ ]:

    
import numpy as np
import logging
import os

import astropy.units as u
import astropy.constants as const
import fastimgproto.imager as imager
import fastimgproto.visibility as visibility

from astropy.coordinates import Angle, SkyCoord, AltAz, EarthLocation
from astropy.time import Time
from fastimgproto.gridder.conv_funcs import GaussianSinc
from fastimgproto.skymodel.helpers import SkyRegion, SkySource
from fastimgproto.sourcefind.image import SourceFindImage
from fastimgproto.telescope.readymade import Meerkat



In [ ]:

    
import matplotlib.pyplot as plt
import matplotlib.image as mpimg
%matplotlib inline
# Plot image pixels in cartesian ordering (i.e. y-positive == upwards):
plt.rcParams['image.origin'] = 'lower'
# Make plots bigger
plt.rcParams['figure.figsize'] = 10, 10

Set up the telescope



In [ ]:

    
telescope = Meerkat()



In [ ]:

    
pointing_centre = SkyCoord(0 * u.deg, -30 * u.deg)
obs_central_frequency = 3. * u.GHz
wavelength = const.c / obs_central_frequency
transit_time = telescope.next_transit(pointing_centre.ra,
                                      start_time=Time('2017-01-01'))

Choose the observation times / sampling regime



In [ ]:

    
nstep=10
obs_times = transit_time + np.linspace(-1, 1, nstep) * u.hr
print("Generating UVW-baselines for {} timesteps".format(nstep))
uvw_m = telescope.uvw_tracking_skycoord(pointing_centre, obs_times)
# From here on we use UVW as multiples of wavelength, lambda:
uvw_lambda = (uvw_m / wavelength).to(u.dimensionless_unscaled).value

Configure the sources to simulate

Steady sources - one at pointing centre, one at a small offset:



In [ ]:

    
# Additional source to North-East of pointing centre
extra_src_position = SkyCoord(ra=pointing_centre.ra + 0.01 * u.deg,
                              dec=pointing_centre.dec + 0.01 * u.deg, )

steady_sources = [
    SkySource(pointing_centre, flux=1 * u.Jy),
    SkySource(extra_src_position, flux=0.4 * u.Jy),
]

And a transient source



In [ ]:

    
transient_posn = SkyCoord(
    ra=pointing_centre.ra - 0.05 * u.deg,
    dec=pointing_centre.dec - 0.05 * u.deg)
transient_sources = [
    SkySource(position=transient_posn, flux=0.5 * u.Jy),
]



In [ ]:

    
model_sources = steady_sources



In [ ]:

    
data_sources = steady_sources
# data_sources = steady_sources + transient_sources
# Simulate some jitter (positional offsets) due to ionospheric effects (or pointing-errors or whatever)
# (ra, dec)

jitter = (0.5*u.arcsec, -0.5*u.arcsec)
jitter_sources = []
for src in data_sources:
    jitter_src = SkySource(position=SkyCoord(ra = src.position.ra + jitter[0], 
                                             dec= src.position.dec+jitter[1]),
                           flux = src.flux
                          )
    jitter_sources.append(jitter_src)



In [ ]:

    
data_sources



In [ ]:

    
jitter_sources

Simulate visibilities, given the source-list and UVW baseline info:



In [ ]:

    
model_vis = visibility.visibilities_for_source_list(
    pointing_centre,
    source_list = model_sources, 
    uvw = uvw_lambda)



In [ ]:

    
# Simulate incoming data; includes transient sources, noise:
baseline_noise_level = 0.1 * u.Jy

data_vis = visibility.visibilities_for_source_list(
    pointing_centre,
#     source_list = data_sources, 
    source_list = jitter_sources,
    uvw = uvw_lambda)

data_vis = visibility.add_gaussian_noise(baseline_noise_level, data_vis)

residual_vis = data_vis - model_vis

Configure and run the imager



In [ ]:

    
image_size=1024 * u.pixel
cell_size=1 * u.arcsecond



In [ ]:

    
kernel_support = 3
kernel_func = GaussianSinc(trunc=kernel_support)



In [ ]:

    
image, beam = imager.image_visibilities(
#     model_vis, 
#     data_vis,
    residual_vis,
    uvw_lambda,
    image_size=image_size,
    cell_size=cell_size,
    kernel_func=kernel_func,
    kernel_support=kernel_support,
    kernel_exact=True,
    kernel_oversampling=None
    )
image = np.real(image)
beam = np.real(beam)

Plot the image and beam-pattern



In [ ]:

    
# %matplotlib notebook
%matplotlib inline
fig, axes = plt.subplots(ncols=2, figsize=(12,8))
clim = (-0.1, 0.7)
# xlim = (250,750)
xlim = (450,550)
# xlim = (550,800)
ylim = xlim

img_ax, bm_ax = axes
im_plot = img_ax.imshow(image,clim=clim)
bm_ax.imshow(beam,clim=clim)

img_ax.set_xlim(*xlim)
img_ax.set_ylim(*ylim)
img_ax.set_title('image')

x_range = xlim[1]-xlim[0]
y_range = ylim[1]-ylim[0]
beam_xlim = ( beam.shape[1]/2 - x_range/2, beam.shape[1]/2 + x_range/2)
beam_ylim = ( beam.shape[0]/2 - y_range/2, beam.shape[0]/2 + y_range/2)
bm_ax.set_xlim(beam_xlim)
bm_ax.set_ylim(beam_ylim)
bm_ax.set_title('beam (matched zoom)')

fig.subplots_adjust(right=0.8)
cbar_ax = fig.add_axes([0.85, 0.15, 0.05, 0.7])
fig.colorbar(im_plot, cax=cbar_ax)
beam_xlim, beam_ylim

Export the image-data to a FITS-file for closer inspection:



In [ ]:

    
import astropy.io.fits as fits
hdu = fits.PrimaryHDU(image)
# hdu.writeto('image.fits')

Sourcefinding



In [ ]:

    
from fastimgproto.sourcefind.image import SourceFindImage



In [ ]:

    
detection_n_sigma=30
analysis_n_sigma=15
sfimage = SourceFindImage(data=np.real(image),
                          detection_n_sigma=detection_n_sigma,
                          analysis_n_sigma=analysis_n_sigma,
                          )



In [ ]:

    
sfimage.islands



In [ ]:

    
src = sfimage.islands[0]

Plot a zoom-in on the first detected source:



In [ ]:

    
fig, ax1 = plt.subplots(1,1)
ax1.imshow(image)
half_width = 75
xlims = int(src.xbar)-half_width, int(src.xbar)+half_width
ylims = int(src.ybar)-half_width, int(src.ybar)+half_width
ax1.set_xlim(xlims)
ax1.set_ylim(ylims)
for src in sfimage.islands:
    ax1.axvline(src.xbar, ls=':')
    ax1.axhline(src.ybar, ls=':')