/
ftplane_convolution.py
116 lines (94 loc) · 3.44 KB
/
ftplane_convolution.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
from spectral_cube import SpectralCube, VaryingResolutionSpectralCube
from radio_beam import Beam
import numpy as np
import astropy.units as u
import astropy.utils.console as console
import copy
import warnings
def ftconvolve(ImageIn, major = 1.0, minor = 1.0,
angle = 0.0):
NanMaskFlag = False
nanmask = np.isnan(ImageIn)
image = np.copy(ImageIn)
if np.any(nanmask):
wtimg = np.ones_like(ImageIn)
NanMaskFlag = True
image[nanmask] = 0.0
wtimg[nanmask] = 0.0
ftwtimg = np.fft.fftn(wtimg)
ftimg = np.fft.fftn(image)
if major == 0.0:
sigmau = np.inf
else:
sigmau = 0.5/(np.pi * major)
if minor == 0.0:
sigmav = np.inf
else:
sigmav = 0.5/(np.pi * minor)
FTPA = angle + np.pi
a = 0.5 * (np.cos(FTPA)**2 / sigmau**2 +
np.sin(FTPA)**2 / sigmav**2)
c = 0.5 * (np.sin(FTPA)**2 / sigmau**2 +
np.cos(FTPA)**2 / sigmav**2)
b = 0.25 * np.sin(2 * FTPA) * (1.0 / sigmav**2 - 1.0 / sigmau**2)
vv, uu = np.meshgrid(np.fft.fftfreq(ftimg.shape[0]),
np.fft.fftfreq(ftimg.shape[1]),
indexing='ij')
FTkernel = np.exp(-a*(uu)**2 -c*(vv)**2 +2*b*(uu*vv))
ConvolvedImage = (np.fft.ifftn(ftimg * FTkernel)).real
if NanMaskFlag:
ConvolvedMask = (np.fft.ifftn(ftwtimg * FTkernel)).real
ConvolvedImage /= ConvolvedMask
ConvolvedImage[nanmask] = np.nan
return(ConvolvedImage)
def MakeRoundBeam(incube, outfile=None):
'''
This takes a FITS file or a SpectralCube and outputs a cube with a
single round beam of maximum size. No provision for robust
statistics is currently made so a single plane with a very large
beam will end up being the target.
Parameters
----------
filename : `string` or `SpectralCube`
Input spectral cube
Returns
-------
cube : `SpectralCube`
'''
if isinstance(incube,str):
cube = SpectralCube.read(incube)
if isinstance(incube, VaryingResolutionSpectralCube):
cube = incube
if not isinstance(cube, VaryingResolutionSpectralCube):
warnings.warn("No information about multiple beams")
return(None)
beams = cube.beams
major = np.array([bm.major.to(u.deg).value for bm in beams])
minor = np.array([bm.minor.to(u.deg).value for bm in beams])
PAs = np.array([bm.pa.to(u.rad).value for bm in beams])
target = np.array(major.max())
# Let's assume square pixels
pixsize = cube.wcs.pixel_scale_matrix[1,1]
output = np.zeros(cube.shape)
fwhm2sigma = np.sqrt(8*np.log(2))
with console.ProgressBar(cube.shape[0]) as bar:
for ii, plane in enumerate(cube.filled_data[:]):
majpix = (target**2 - major[ii]**2)**0.5 /\
pixsize / fwhm2sigma
minpix = (target**2 - minor[ii]**2)**0.5 /\
pixsize / fwhm2sigma
output[ii,:,:] = ftconvolve(plane,
major = majpix,
minor = minpix,
angle = PAs[ii])
bar.update()
hdr = copy.copy(cube.header)
hdr['CASAMBM'] = False
hdr['BMAJ'] = float(target)
hdr['BMIN'] = float(target)
hdr['BPA'] = 0.0
outcube = SpectralCube(output, cube.wcs, header=hdr)
if outfile:
outcube.write(outfile)
return None
return(outcube)