-
Notifications
You must be signed in to change notification settings - Fork 10
/
fits2sky.py
executable file
·204 lines (175 loc) · 7.21 KB
/
fits2sky.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
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
#!/usr/bin/env python
# -*- coding: utf-8 -*-
#
# Copyright (C) 2019 - Francesco de Gasperin
#
# This program is free software; you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
# the Free Software Foundation; either version 2 of the License, or
# (at your option) any later version.
#
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
"""
Script to make a sky model from fits model images
"""
import argparse
from argparse import RawTextHelpFormatter
from astropy.io import fits
from astropy import wcs
import numpy as np
import scipy.interpolate
import casacore.tables as pt
import sys
import os
import glob
def ra2hhmmss(deg):
"""Convert RA coordinate (in degrees) to HH MM SS"""
from math import modf
if deg < 0:
deg += 360.0
x, hh = modf(deg/15.)
x, mm = modf(x*60)
ss = x*60
return (int(hh), int(mm), ss)
def dec2ddmmss(deg):
"""Convert DEC coordinate (in degrees) to DD MM SS"""
from math import modf
sign = (-1 if deg < 0 else 1)
x, dd = modf(abs(deg))
x, ma = modf(x*60)
sa = x*60
return (int(dd), int(ma), sa, sign)
def convert_radec_str(ra, dec):
"""Takes ra, dec in degrees and returns makesourcedb strings"""
ra = ra2hhmmss(ra)
sra = str(ra[0]).zfill(2)+':'+str(ra[1]).zfill(2)+':'+str("%.3f" % (ra[2])).zfill(6)
dec = dec2ddmmss(dec)
decsign = ('-' if dec[3] < 0 else '+')
sdec = decsign+str(dec[0]).zfill(2)+'.'+str(dec[1]).zfill(2)+'.'+str("%.3f" % (dec[2])).zfill(6)
return sra, sdec
def main(fits_model_root, skymodel, ref_freq='60e6', fits_mask=None, min_peak_flux_jy=0.001,
max_residual_jy=0.00, interp='linear'):
"""
Make a makesourcedb sky model for input MS from WSClean fits model images
Parameters
----------
fits_model_root : str
Root name of WSClean fits model files (without the "-XXXX-model.fits" part)
skymodel : str
Filename of the output makesourcedb sky model
ref_freq : float, optional
Reference freq of the output catalogue in Hz
fits_mask : str, optional
Filename of fits mask
min_peak_flux_jy : float, optional
Minimum absolute value of flux in Jy of a source in lowest-frequency model image
to include in output model
max_residual_jy : float, optional
Maximum acceptible total residual absolute flux in Jy
interp : str, optional
Interpolation method. Can be any supported by scipy.interpolate.interp1d:
'linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic'
"""
min_peak_flux_jy = float(min_peak_flux_jy)
max_residual_jy = float(max_residual_jy)
if type(fits_mask) is str:
if fits_mask.lower() == 'none':
fits_mask = None
# Find model images: look first for channel images and MFS image
fits_models = glob.glob(fits_model_root+'-*-model.fits')
if len(fits_models) > 0:
# Get the MFS image
mfs_model = fits_model_root+'-MFS-model.fits'
else:
# No channels images found, so look for non-MFS images
fits_models = glob.glob(fits_model_root+'-model.fits')
mfs_model = None
if len(fits_models) == 0:
print('ERROR: no model images found')
sys.exit(1)
# Read in model images
freqs = []
model_images = []
for f in fits_models:
# Get the frequency info
hdr = fits.getheader(f, 0, ignore_missing_end=True)
freqs.append(hdr['CRVAL3']) # Hz
model_images.append(fits.getdata(f, 0, ignore_missing_end=True))
w = wcs.WCS(hdr)
# Read in MFS image
if mfs_model is None:
mfs_model = fits_models[0]
mfs_image = fits.getdata(mfs_model, 0, ignore_missing_end=True)
# Sort by freq
sorted_ind = np.argsort(freqs)
freqs = np.array(freqs)[sorted_ind]
fits_models = np.array(fits_models)[sorted_ind]
model_images = np.array(model_images)[sorted_ind]
# Find pixels that meet the flux cut (and are in the mask, if given)
if fits_mask is not None:
if fits_mask.lower() == 'empty':
# Handle case in which no sources were found during masking
nonzero_ind = [[], []]
else:
mask = fits.getdata(fits_mask, 0, ignore_missing_end=True)
nonzero_ind = np.where((np.abs(mfs_image) > min_peak_flux_jy) & (mask > 0))
else:
nonzero_ind = np.where(np.abs(mfs_image) > min_peak_flux_jy)
# Interpolate the fluxes to the frequency of the MS
nsources = len(nonzero_ind[0])
fluxes = []
names = []
ras = []
decs = []
for i in range(nsources):
index = [nonzero_ind[j][i] for j in range(4)]
index.reverse() # change to WCS coords
ras.append(w.wcs_pix2world(np.array([index]), 0, ra_dec_order=True)[0][0])
decs.append(w.wcs_pix2world(np.array([index]), 0, ra_dec_order=True)[0][1])
names.append('cc{}'.format(i))
index.reverse() # change back to image coords
flux_array = np.array([im[tuple(index)] for im in model_images])
# If MS frequency lies outside range, just use nearest freq
if ref_freq < freqs[0]:
flux = flux_array[0]
elif ref_freq > freqs[-1]:
flux = flux_array[-1]
else:
# Otherwise interpolate
flux = scipy.interpolate.interp1d(freqs, flux_array, kind=interp)(ref_freq)
fluxes.append(flux)
# Remove sources until we reach the desired residual
if len(fluxes) > 0:
total_flux = np.sum(np.abs(fluxes))
keep_ind = np.where(np.abs(fluxes) > min_peak_flux_jy)
while (total_flux - np.sum(np.abs(np.array(fluxes)[keep_ind]))) < max_residual_jy:
min_peak_flux_jy *= 1.1
keep_ind = np.where(np.abs(fluxes) > min_peak_flux_jy)
if len(keep_ind[0]) < 50:
# keep up to 50 sources regardless of the residual
break
fluxes = np.array(fluxes)[keep_ind]
ras = np.array(ras)[keep_ind]
decs = np.array(decs)[keep_ind]
names = np.array(names)[keep_ind]
# Write sky model
with open(skymodel, 'w') as outfile:
outfile.write('FORMAT = Name, Type, Ra, Dec, I, Q, U, V, ReferenceFrequency\n')
for name, ra, dec, flux in zip(names, ras, decs, fluxes):
ra_str, dec_str = convert_radec_str(ra, dec)
outfile.write('{0}, POINT, {1}, {2}, {3}, 0.0, 0.0, 0.0, {4}\n'
.format(name, ra_str, dec_str, flux, ref_freq))
if __name__ == '__main__':
descriptiontext = "Make a makesourcedb sky model from WSClean fits model images.\n"
parser = argparse.ArgumentParser(description=descriptiontext, formatter_class=RawTextHelpFormatter)
parser.add_argument('fits_model_root', help='Root of model images')
parser.add_argument('skymodel', help='Filename of output sky model')
args = parser.parse_args()
main(args.fits_model_root, args.skymodel)