Feathering

This talk can be found at http://small.cat/hat until next week. Please follow along on your laptop or phone.

Goals:

  • Describe and demonstrate feathering process
  • Show some of the free parameters that can be selected when feathering
  • Show various visualizations that can be used to assess goodness-of-merge
    • PSD plots
    • Residuals from theoretical
  • Demonstrate limitations of image combination
  • Show tools for simulating realistic synthetic data

What is feathering?

Fourier-space (UV-space) weighted averaging of two data sets.

The weight function is selected to match the Fourier transform of the single-dish telescope's beam.

Benefits & Costs of feathering

The good:

  • Simple
  • Applicable to non-interferometer (e.g., bolometer + other) data

The bad:

  • Deceptively simple
  • Ignores spectral information (see SDINT, TP2VIS talks)
  • Does not fix interferometric image reconstruction problems
  • Tricky to balance weights for correctness vs optimal noise

radio-astro-tools

http://radio-astro-tools.github.io

Includes tools for cubes (https://spectral-cube.readthedocs.io), beam metadata handling (https://github.com/radio-astro-tools/radio-beam), and feathering (https://github.com/radio-astro-tools/uvcombine/)

Feathering (combination) tests

This notebook presents a series of experiments in single-dish + interferometer combination on "realistic" data.

We're "observing" at 2mm, so a 12m dish has a FWHM=40" Gaussian beam and a 9m baseline has a sharp cutoff at 56"

This presentation deals mostly with idealized cases; there is no primary beam degradation included in the simulations.

Requirements for this work: turbustat generates our synthetic data and helps with power-spectral-density (PSD) plotting. astropy.convolution provides access to convolution tools, and uvcombine is our python-only implementation of feather.

https://turbustat.readthedocs.io/en/latest/, especially https://turbustat.readthedocs.io/en/latest/generating_test_data.html

from turbustat.simulator.gen_field import make_extended
from turbustat.statistics import psds
from astropy import convolution, units as u
import numpy as np
from uvcombine.uvcombine import feather_kernel, fftmerge

We create a synthetic power-law power-spectrum image. This sort of image is typical of a dust image of the Galactic plane, for example.

# create an input image with specified parameters
# (this can later be modified - it will be good to examine the effects of
# different power laws, different types of input...)
# We're assuming a scale of 1"/pixel for this example
imsize = 512
powerlaw = 3
im = make_extended(imsize=imsize, powerlaw=powerlaw, randomseed=0)
# the real sky is positive, so we subtract the minimum to force the overall image positive
im = im - im.min()

Input Image Visualization

This is the input image along with its histogram.

The power spectrum of the input image (set to be $\alpha=3$), verifying that the turbustat code works

Next, we create our simulated interferometer by creating a UV domain and selecting which pixels in that domain will be part of our telescope. This process creates an ideal interferometer.

# set up the grid
ygrid, xgrid = np.indices(im.shape, dtype='float')
rr = ((xgrid-im.shape[1]/2)**2+(ygrid-im.shape[0]/2)**2)**0.5

# Create a UV sampling mask.
# This removes all large-angular scale (r<8) features *in UV space* and all
# small angular scales.
# In fourier space, r=0 corresponds to the DC component
# r=1 corresponds to the full map (one period over that map)
# r=256 is the smallest angular scale, which is 2 pixels
# We're assuming a pixel scale of 1" / pixel
# therefore 56" corresponds to 9m at 2mm (i.e., nearly the closest spacing possible for 7m)
# We cut off the "interferometer" at 2.5" resolution
largest_scale = 56.*u.arcsec
smallest_scale = 2.5*u.arcsec
pixscale = 1*u.arcsec
image_scale = im.shape[0]*pixscale # assume symmetric (default=256)
ring = (rr>=(image_scale/largest_scale)) & (rr<=(image_scale/smallest_scale))

Synthetic Perfect Interferometer

The synthetic interferometer's UV coverage map (it's a perfect interferometer)

Next, we create the interferometric map by multiplying our interferometer mask by the fourier transform of the sky data

# create the interferometric map by removing both large and small angular
# scales in fourier space
imfft = np.fft.fft2(im)
imfft_interferometered = imfft * np.fft.fftshift(ring)
im_interferometered = np.fft.ifft2(imfft_interferometered)

The interferometric image does not preserve total flux, as expected. Note that the mean of the histogram is shifted.

The residual of the original image minus the interferometrically observed image. The large scales and noise are preserved.