Feathering¶
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.
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"
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.
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
# im = make_extended(imsize=imsize, powerlaw=1.5, randomseed=0)
# # the real sky is positive, so we subtract the minimum to force the overall image positive
# im = im - im.min()
# get data
import requests
from astropy.io import fits
if not os.path.exists('orionA-C-250.fits.gz'):
req = requests.get("http://www.herschel.fr/cea/gouldbelt/en/archives/orionA/orionA-C-250.fits.gz", stream=True)
with open("orionA-C-250.fits.gz", "wb") as fh:
fh.write(req.content)
fh = fits.open('orionA-C-250.fits.gz')
im = fh[0].data[512:1024, 512:1024]
Input Image Visualization¶
This is the input image along with its histogram.
The power spectrum of the input image (still pretty close to $\alpha=1.5$)
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.
Synthetic Single Dish¶
The single dish map is just a convolution of the original data with a Gaussian beam. It preserves flux but loses small scales.
# create the single-dish map by convolving the image with a FWHM=40" kernel
# (this interpretation is much easier than the sharp-edged stuff in fourier
# space because the kernel is created in real space)
lowresfwhm = 40*u.arcsec
singledish_kernel = convolution.Gaussian2DKernel(lowresfwhm/pixscale/2.35, x_size=im.shape[1], y_size=im.shape[0])
singledish_kernel_fft = np.fft.fft2(singledish_kernel)
singledish_im = convolution.convolve_fft(im,
kernel=singledish_kernel,
boundary='fill',
fill_value=im.mean())
The single-dish image and its histogram
The single dish in Fourier space¶
We show the single dish beam in Fourier space with the interferometer coverage range overlaid
Feathering¶
Feathering is the combination of the single-dish image with the interferometric image in the UV domain.
In the uvcombine package, this is handled by uvcombine.feather_simple. However, we show the components of that function here.
For comparison, CASA's feather takes these inputs:
# feather :: Combine two images using their Fourier transforms
imagename = '' # Name of output feathered image
highres = '' # Name of high resolution (interferometer) image
lowres = '' # Name of low resolution (single dish) image
sdfactor = 1.0 # Scale factor to apply to Single Dish image
effdishdiam = -1.0 # New effective SingleDish diameter to use in m
lowpassfiltersd = False # Filter out the high spatial frequencies of the SD imageFirst, we define the FWHM of the low-resolution (single-dish) image, which defines the effective cutoff point between the interferometer and the single dish. lowresfwhm is equivalent to effdishdiam from CASA, albeit in different units.
The kernels weight arrays for the single-dish and interferometer data. They are the fourier transforms of the low-resolution beam and (1-that kernel), respectively.
# pixel scale can be interpreted as "arcseconds"
# then, fwhm=40 means a beam fwhm of 40"
pixscale = 1*u.arcsec
lowresfwhm = 40*u.arcsec
nax1,nax2 = im.shape
kfft, ikfft = feather_kernel(nax2, nax1, lowresfwhm, pixscale,)
Then we specify a few parameters that are not all available in CASA. CASA's lowpassfiltersd is equivalent to our replace_hires, and their sdfactor is our lowresscalefactor. The other parameters, highresscalefactor and deconvSD are unavailable in CASA.
# Feather the interferometer and single dish data back together
# This uses the naive assumptions that CASA uses
# However, there are a few flags that can be played with.
# None of them do a whole lot, though there are good theoretical
# reasons to attempt them.
im_hi = im_interferometered.real
im_low = singledish_im
lowresscalefactor=1
replace_hires = False
lowpassfilterSD = False
deconvSD = False
highresscalefactor=1
fftsum, combo = fftmerge(kfft, ikfft, im_hi*highresscalefactor,
im_low*lowresscalefactor,
replace_hires=replace_hires,
lowpassfilterSD=lowpassfilterSD,
deconvSD=deconvSD,
)
The feathered data set looks pretty good.
This image looks pretty close to the original, but the peaks and valleys certainly are not recovered (the contrast is reduced compared to the original).
We then compare the feathered data to the input image.
The difference between the input image and the feathered image shows the remainder artifacts.
If we repeat the same experiment again, but with the shortest baseline at 12m instead of 9m, the results are noticeably worse:
largest_scale = 42.*u.arcsec # (1.22 * 2*u.mm/(12*u.m)).to(u.arcsec, u.dimensionless_angles())
smallest_scale = 2.5*u.arcsec
image_scale = im.shape[0]*pixscale # assume symmetric (default=256)
ring = (rr>=(image_scale/largest_scale)) & (rr<=(image_scale/smallest_scale))
# 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)
im_hi = im_interferometered.real
lowresscalefactor=1
replace_hires = False
lowpassfilterSD = False
deconvSD = False
highresscalefactor=1
fftsum, combo = fftmerge(kfft, ikfft, im_hi*highresscalefactor,
im_low*lowresscalefactor,
replace_hires=replace_hires,
lowpassfilterSD=lowpassfilterSD,
deconvSD=deconvSD,
)
The feathered image with 12m shortest baselines instead of 9m
The side-by-side images comparison again, but with 12m instead of 9m shortest baselines
The difference image (original - feathered w/12m shortest baselines)
It is more helpful to look at the difference in power spectra. Note that the axes are log-scaled.
The components that go in to the feather can help clarify the picture
Different combinations of parameters can yield very different results.
We have 8 parameter combinations:
- Replace singledish-interferometer overlap w/interferometer, or average them
- Filter out the small angular scales from the single dish or don't
- Deconvolve the single dish (direct deconvolution) or don't
The next slides are a walkthrough of parameter exploration with different effective dish sizes.
First, we have a 12m single-dish combined with an interferometer with shortest baseline length of 12m.
import imp
import uvcombine.plot_utilities
imp.reload(uvcombine.plot_utilities)
from uvcombine.plot_utilities import compare_parameters_feather_simple
Residuals are in blue.
A short-spacing limit of 12m for the interferometer and dish size of 12m is a bad case, but not the worst.
ALMA's shortest baselines are 14.6m (main array) and 8.7m (7m array). L05, the 5th percentile baseline length, is 21.4m (9.1m) for the 12m (7m) array.
A realistic case for ALMA is then a 9m shortest baseline and a 12m effective single dish.
Residuals are in blue.
This realistic case is still bad. The "best case", in which we have good UV overlap (e.g., a 24m dish instead of a 12m), actually can get good theoretical recovery
Residuals are in blue.
CASA's default parameters are fine, but for best performance, the single-dish image should be deconvolved.
The appropriate choice of feather parameters depends (mildly) on the UV overlap:
- If the single dish is substantially bigger than the shortest baseline, CASA's defaults or replacing short-spacing with deconvolved single-dish both work well
- If the single dish is comparable to the shortest baseline, the best results come from deconvolving the single-dish data and weighted-averaging them with the interferometric
These feather experiments represent the absolute best-case scenario. They should be used as references for comparison with any other combination technique.
Non-ideal cases¶
Besides simple UV coverage problems (which are expensive to fix and usually out of our control), there are other issues:
- Relative calibration
- Relative sensitivity
Relative calibration is an issue if the data simply aren't calibrated perfectly (systematic uncertainties are usually 5-15%) or if the single-dish data come from a different frequency range (e.g., using wide-band bolometer data to fill in the continuum zero-spacing).
We can simulate calibration error using the highresscalefactor and lowresscalefactor parameters, which are included to correct for these errors
lowresscalefactor = 1.50
highresscalefactor = 1.0
fftsum, combo = fftmerge(kfft, ikfft, im_hi*highresscalefactor,
im_low*lowresscalefactor,
replace_hires=replace_hires,
lowpassfilterSD=lowpassfilterSD,
deconvSD=deconvSD,
)
lowresscalefactor = 1.0
highresscalefactor = 1.50
fftsum, combo = fftmerge(kfft, ikfft, im_hi*highresscalefactor,
im_low*lowresscalefactor,
replace_hires=replace_hires,
lowpassfilterSD=lowpassfilterSD,
deconvSD=deconvSD,
)
Relative sensitivity¶
The cases we've treated above assume that the interferometer and single-dish telescope have infinite sensitivity; the "noise" in the image is actually part of the sky. So, what happens if we add observational (as opposed to astrophysical) noise?
import radio_beam
import astropy.utils.misc
lowresscalefactor = 1.0
highresscalefactor = 1.0
sd_beam = radio_beam.Beam(lowresfwhm)
sd_beam_volume = (sd_beam.sr / pixscale**2).decompose()
with astropy.utils.misc.NumpyRNGContext(0):
noise_sd = convolution.convolve_fft(singledish_im.max() * 0.0002 * np.random.randn(*singledish_im.shape) * sd_beam_volume,
sd_beam.as_kernel(pixscale))
noisy_sd = (singledish_im + noise_sd)
fftsum, combo = fftmerge(kfft, ikfft, im_hi*highresscalefactor,
noisy_sd*lowresscalefactor,
replace_hires=replace_hires,
lowpassfilterSD=lowpassfilterSD,
deconvSD=deconvSD,
)
Now we'll do the same experiment with noise added to the interferometric image too
lowresscalefactor = 1.0
highresscalefactor = 1.0
intf_beam = radio_beam.Beam(u.Quantity(smallest_scale, u.arcsec))
intf_beam_volume = (intf_beam.sr / pixscale**2).decompose()
assert intf_beam_volume < 100
with astropy.utils.misc.NumpyRNGContext(0):
noise_intf = convolution.convolve_fft(im_interferometered.real.max() * 0.02 * np.random.randn(*im_interferometered.shape) * intf_beam_volume,
intf_beam.as_kernel(pixscale))
noisy_intf = (im_interferometered.real + noise_intf)
fftsum, combo = fftmerge(kfft, ikfft, noisy_intf*highresscalefactor,
noisy_sd*lowresscalefactor,
replace_hires=replace_hires,
lowpassfilterSD=lowpassfilterSD,
deconvSD=deconvSD,
)