Commits
Pam Harris authored 300af8941f5 Merge
1 + | ########################################################################## |
2 + | # test_stk_vla_users_continuum_from_SDM.py |
3 + | # |
4 + | # Copyright (C) 2018 |
5 + | # Associated Universities, Inc. Washington DC, USA. |
6 + | # |
7 + | # This script is free software; you can redistribute it and/or modify it |
8 + | # under the terms of the GNU Library General Public License as published by |
9 + | # the Free Software Foundation; either version 2 of the License, or (at your |
10 + | # option) any later version. |
11 + | # |
12 + | # This library is distributed in the hope that it will be useful, but WITHOUT |
13 + | # ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
14 + | # FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public |
15 + | # License for more details. |
16 + | # |
17 + | # https://open-jira.nrao.edu/browse/CASR-662 |
18 + | # https://open-jira.nrao.edu/browse/CAS-14036 |
19 + | # |
20 + | # |
21 + | ########################################################################## |
22 + | |
23 + | ''' |
24 + | VLA Users Stakeholder Continuum Test |
25 + | |
26 + | Based on 3C 391 CASAguide Continuum Tutorial for CASA Version 6.2.0 |
27 + | https://casaguides.nrao.edu/index.php?title=VLA_Continuum_Tutorial_3C391-CASA6.2.0 |
28 + | |
29 + | Test list |
30 + | - flag fraction and MS nrows after import, split and editing steps |
31 + | - antenna position corrections retrieved by gencal |
32 + | - fluxscale input table and bootstrapping results |
33 + | - imstat values from the final pbcor image |
34 + | - imstat values from the final residual image |
35 + | ''' |
36 + | |
37 + | ########################################################################## |
38 + | ########################################################################## |
39 + | |
40 + | # Imports # |
41 + | import os |
42 + | import unittest |
43 + | import numpy as np |
44 + | import shutil |
45 + | from glob import glob |
46 + | |
47 + | from casatestutils import stats_dict |
48 + | from casatestutils import generate_weblog |
49 + | from casatestutils import add_to_dict |
50 + | |
51 + | from casatools import ctsys |
52 + | from casatools import table as tbtool |
53 + | from casatasks import casalog, importasdm, listobs, flagdata, gencal, setjy, gaincal, bandpass |
54 + | from casatasks import fluxscale, applycal, split, statwt, tclean, impbcor, imstat |
55 + | from casatasks.private.parallel.parallel_task_helper import ParallelTaskHelper |
56 + | |
57 + | # location of data |
58 + | data_path = ctsys.resolve('stakeholder/vla/') |
59 + | |
60 + | clean_mask = "#CRTFv0 CASA Region Text Format version 0\npoly [[18:49:21.55474, -000.52.32.9190], [18:49:29.15960, -000.53.10.9387], [18:49:31.21928, -000.53.48.9583], [18:49:31.53618, -000.54.38.8592], [18:49:35.33866, -000.55.02.6210], [18:49:40.09177, -000.55.21.6298], [18:49:40.72560, -000.56.25.7880], [18:49:41.04265, -000.58.36.4810], [18:49:38.82455, -000.59.14.5014], [18:49:31.37791, -000.59.40.6415], [18:49:23.29749, -000.59.21.6320], [18:49:24.08969, -000.58.17.4736], [18:49:23.29750, -000.57.56.0874], [18:49:15.37558, -000.57.53.7105], [18:49:11.73154, -000.57.13.3138], [18:49:11.41470, -000.56.37.6701], [18:49:16.16792, -000.54.07.9680], [18:49:19.65353, -000.52.42.4239]] coord=J2000, corr=[I]" |
61 + | |
62 + | |
63 + | ############################################## |
64 + | test_dict = {} |
65 + | #### Tests #### |
66 + | class Test_vla_users_continuum(unittest.TestCase): |
67 + | |
68 + | def setUp(self): |
69 + | self._mytb = tbtool() |
70 + | self.refversion='6.2.0.124' |
71 + | self.sdmfile = 'TDEM0001_sb1218006_1.55310.33439732639' |
72 + | os.symlink(data_path+self.sdmfile, self.sdmfile) |
73 + | self.maskfile = '3c391_clean_mask.crtf' |
74 + | # os.symlink(data_path+self.maskfile, self.maskfile) |
75 + | self.writeMask() |
76 + | # Control if running in parallel with mpicasa |
77 + | self.parallel = False |
78 + | # Control if creating Multi-MS with mpicasa |
79 + | self.mms = False |
80 + | if ParallelTaskHelper.isMPIEnabled(): |
81 + | casalog.post('MPI is enabled. Will run in parallel using Multi-MS for calibration and tclean') |
82 + | self.parallel = True |
83 + | self.mms = True |
84 + | |
85 + | def tearDown(self): |
86 + | generate_weblog("vla_users_continuum_from_SDM",test_dict) |
87 + | self._mytb.done() |
88 + | self.delData() |
89 + | |
90 + | def writeMask(self): |
91 + | with open(self.maskfile,'w') as out1: |
92 + | out1.write(clean_mask) |
93 + | |
94 + | def delData(self): |
95 + | os.unlink( self.sdmfile ) |
96 + | os.unlink( self.maskfile ) |
97 + | os.system('rm -rf '+self.sdmfile+'*') |
98 + | del_files = glob('3c391_ctm*') |
99 + | for f in del_files: |
100 + | shutil.rmtree(f) |
101 + | |
102 + | test_dict) | (
103 + | def test_continuum_from_SDM(self): |
104 + | """VLA Stakeholders tests: Based on 3C 391 CASAguide Continuum Tutorial""" |
105 + | sdmname = self.sdmfile |
106 + | msname = '3c391_ctm_mosaic_10s_spw0.ms' |
107 + | msname_split = '3c391_ctm_mosaic_spw0.ms' |
108 + | |
109 + | test_name = self._testMethodName |
110 + | report=[] |
111 + | |
112 + | ## Data Import |
113 + | importasdm(asdm=sdmname, |
114 + | vis=sdmname+'.ms', createmms=self.mms, |
115 + | ocorr_mode='co', asis='Receiver CalAtmosphere', |
116 + | process_caldevice=True, process_pointing=True, savecmds=True, |
117 + | outfile=sdmname+'.flagonline.txt', |
118 + | overwrite=False, bdfflags=False, with_pointing_correction=True) |
119 + | |
120 + | add_to_dict(self, output = test_dict, dataset = sdmname+'.ms') |
121 + | |
122 + | flagdata(vis=sdmname+'.ms', mode='list', |
123 + | inpfile=sdmname+'.flagonline.txt', |
124 + | tbuff=0.0, action='apply', flagbackup=False, savepars=False) |
125 + | |
126 + | flagdata(vis=sdmname+'.ms', mode='clip', clipzeros=True, |
127 + | action='apply', flagbackup=False, savepars=False) |
128 + | |
129 + | flagdata(vis=sdmname+'.ms', mode='shadow', tolerance=0.0, |
130 + | action='apply', flagbackup=False, savepars=False) |
131 + | |
132 + | split(vis=sdmname+'.ms', outputvis=msname, keepflags=False, |
133 + | spw='0', timebin='10s', datacolumn='DATA') |
134 + | |
135 + | add_to_dict(self, output = test_dict, dataset = msname) |
136 + | |
137 + | |
138 + | ## Data Editing |
139 + | flagdata(vis=msname, flagbackup=True, mode='manual', scan='1') |
140 + | |
141 + | flagdata(vis=msname, flagbackup=True, mode='manual', antenna='ea13,ea15') |
142 + | |
143 + | flagdata(vis=msname, mode='quack', quackinterval=10.0, quackmode='beg') |
144 + | |
145 + | flagsummary = flagdata(vis=msname, mode='summary') |
146 + | result, expected, tolerance = 100* flagsummary['spw']['0']['flagged']/flagsummary['spw']['0']['total'], 20.558870499248552, 1.0 |
147 + | msg = f"Expected flag percentage of {expected}, got {result}, tolerance {tolerance}" |
148 + | print(msg); report.append(msg) |
149 + | self.assertAlmostEqual(result, expected, delta=tolerance , msg=msg ) |
150 + | |
151 + | |
152 + | listobs_output = listobs(vis=msname) |
153 + | result, expected = listobs_output['numrecords'], 850358 #845379 |
154 + | msg = f"Expected {expected} rows in split MS, got {result}, tolerance exact" |
155 + | print(msg); report.append(msg) |
156 + | self.assertTrue(result==expected, msg=msg ) |
157 + | |
158 + | |
159 + | ## Calibration |
160 + | gencal(vis=msname,caltable='3c391_ctm_mosaic_10s_spw0.antpos',caltype='antpos') |
161 + | |
162 + | # TEST: compare contents of gencal table |
163 + | self._mytb.open('3c391_ctm_mosaic_10s_spw0.antpos') |
164 + | ant_array=self._mytb.getcol('FPARAM').squeeze().T |
165 + | self._mytb.close() |
166 + | nonzero_ant_array=ant_array[np.logical_not(np.all( ant_array==0, axis=1 ))] |
167 + | |
168 + | result, expected = nonzero_ant_array.shape[0], 14 |
169 + | msg = f"Expected {expected} nonzero antenna position corrections, got {result}, tolerance exact" |
170 + | print(msg); report.append(msg) |
171 + | self.assertTrue(result==expected, msg=msg ) |
172 + | |
173 + | result, expected, tolerance = np.min(nonzero_ant_array), -0.0257, 0.00001 |
174 + | msg = f"Expected minimum antenna position correction of {expected}, got {result}, tolerance {tolerance}" |
175 + | print(msg); report.append(msg) |
176 + | self.assertAlmostEqual(result, expected, delta=tolerance , msg=msg ) |
177 + | |
178 + | result, expected, tolerance = np.max(nonzero_ant_array), 0.0045, 0.00001 |
179 + | msg = f"Expected maximum antenna position correction of {expected}, got {result}, tolerance {tolerance}" |
180 + | print(msg); report.append(msg) |
181 + | self.assertAlmostEqual(result, expected, delta=tolerance , msg=msg ) |
182 + | |
183 + | |
184 + | setjy(vis=msname,field='J1331+3030',standard='Perley-Butler 2017', |
185 + | model='3C286_C.im',usescratch=True,scalebychan=True,spw='') |
186 + | |
187 + | flagdata(vis=msname, |
188 + | flagbackup=True, mode='manual', antenna='ea05') |
189 + | |
190 + | gaincal(vis=msname, caltable='3c391_ctm_mosaic_10s_spw0.G0', |
191 + | field='J1331+3030', refant='ea21', spw='0:27~36', calmode='p', solint='int', |
192 + | minsnr=5, gaintable=['3c391_ctm_mosaic_10s_spw0.antpos']) |
193 + | |
194 + | gaincal(vis=msname,caltable='3c391_ctm_mosaic_10s_spw0.K0', |
195 + | field='J1331+3030',refant='ea21',spw='0:5~58',gaintype='K', |
196 + | solint='inf',combine='scan',minsnr=5, |
197 + | gaintable=['3c391_ctm_mosaic_10s_spw0.antpos', |
198 + | '3c391_ctm_mosaic_10s_spw0.G0']) |
199 + | |
200 + | bandpass(vis=msname,caltable='3c391_ctm_mosaic_10s_spw0.B0', |
201 + | field='J1331+3030',spw='',refant='ea21',combine='scan', |
202 + | solint='inf',bandtype='B', |
203 + | gaintable=['3c391_ctm_mosaic_10s_spw0.antpos', |
204 + | '3c391_ctm_mosaic_10s_spw0.G0', |
205 + | '3c391_ctm_mosaic_10s_spw0.K0']) |
206 + | |
207 + | gaincal(vis=msname,caltable='3c391_ctm_mosaic_10s_spw0.G1', |
208 + | field='J1331+3030',spw='0:5~58', |
209 + | solint='inf',refant='ea21',gaintype='G',calmode='ap',solnorm=False, |
210 + | gaintable=['3c391_ctm_mosaic_10s_spw0.antpos', |
211 + | '3c391_ctm_mosaic_10s_spw0.K0', |
212 + | '3c391_ctm_mosaic_10s_spw0.B0'], |
213 + | interp=['','','nearest']) |
214 + | |
215 + | gaincal(vis=msname,caltable='3c391_ctm_mosaic_10s_spw0.G1', |
216 + | field='J1822-0938', |
217 + | spw='0:5~58',solint='inf',refant='ea21',gaintype='G',calmode='ap', |
218 + | gaintable=['3c391_ctm_mosaic_10s_spw0.antpos', |
219 + | '3c391_ctm_mosaic_10s_spw0.K0', |
220 + | '3c391_ctm_mosaic_10s_spw0.B0'], |
221 + | append=True) |
222 + | |
223 + | myscale = fluxscale(vis=msname, |
224 + | caltable='3c391_ctm_mosaic_10s_spw0.G1', |
225 + | fluxtable='3c391_ctm_mosaic_10s_spw0.fluxscale1', |
226 + | reference='J1331+3030', |
227 + | transfer=['J1822-0938'], |
228 + | incremental=False) |
229 + | |
230 + | |
231 + | # TEST: compare specific values in return dictionary |
232 + | result, expected = myscale['1']['0']['numSol'][0], 46 |
233 + | msg = f"Expected {expected} input fluxboot solutions, got {result}, tolerance exact" |
234 + | print(msg); report.append(msg) |
235 + | self.assertTrue(result==expected, msg=msg ) |
236 + | |
237 + | result, expected, epsilon = myscale['1']['0']['fluxd'][0], 2.2973563472684653, 0.04 |
238 + | msg = f"Expected J1822-0938 bootstrapped flux of {expected}, got {result}, tolerance {100*epsilon}%" |
239 + | print(msg); report.append(msg) |
240 + | self.assertAlmostEqual(result, expected, delta=expected*epsilon, msg=msg ) |
241 + | |
242 + | |
243 + | applycal(vis=msname, |
244 + | field='J1331+3030', |
245 + | gaintable=['3c391_ctm_mosaic_10s_spw0.antpos', |
246 + | '3c391_ctm_mosaic_10s_spw0.fluxscale1', |
247 + | '3c391_ctm_mosaic_10s_spw0.K0', |
248 + | '3c391_ctm_mosaic_10s_spw0.B0'], |
249 + | gainfield=['','J1331+3030','',''], |
250 + | interp=['','nearest','',''], |
251 + | calwt=False) |
252 + | |
253 + | applycal(vis=msname, |
254 + | field='J1822-0938', |
255 + | gaintable=['3c391_ctm_mosaic_10s_spw0.antpos', |
256 + | '3c391_ctm_mosaic_10s_spw0.fluxscale1', |
257 + | '3c391_ctm_mosaic_10s_spw0.K0', |
258 + | '3c391_ctm_mosaic_10s_spw0.B0'], |
259 + | gainfield=['','J1822-0938','',''], |
260 + | interp=['','nearest','',''], |
261 + | calwt=False) |
262 + | |
263 + | applycal(vis=msname, |
264 + | field='2~8', |
265 + | gaintable=['3c391_ctm_mosaic_10s_spw0.antpos', |
266 + | '3c391_ctm_mosaic_10s_spw0.fluxscale1', |
267 + | '3c391_ctm_mosaic_10s_spw0.K0', |
268 + | '3c391_ctm_mosaic_10s_spw0.B0'], |
269 + | gainfield=['','J1822-0938','',''], |
270 + | interp=['','linear','',''], |
271 + | calwt=False) |
272 + | |
273 + | |
274 + | ## Imaging |
275 + | split(vis=msname,outputvis=msname_split, |
276 + | datacolumn='corrected',field='2~8',correlation='RR,LL') |
277 + | |
278 + | add_to_dict(self, output = test_dict, dataset = msname_split) |
279 + | |
280 + | statwt(vis=msname_split,datacolumn='data') |
281 + | |
282 + | tclean(vis=msname_split,imagename='3c391_ctm_spw0_multiscale', |
283 + | field='',spw='', |
284 + | specmode='mfs', |
285 + | niter=20000, |
286 + | gain=0.1, threshold='1.0mJy', |
287 + | gridder='mosaic', |
288 + | deconvolver='multiscale', |
289 + | scales=[0, 5, 15, 45], smallscalebias=0.9, |
290 + | interactive=False, |
291 + | imsize=[480,480], cell=['2.5arcsec','2.5arcsec'], |
292 + | stokes='I', |
293 + | weighting='briggs',robust=0.5, |
294 + | pbcor=False, |
295 + | mask=self.maskfile, |
296 + | savemodel='none', |
297 + | parallel=self.parallel) |
298 + | |
299 + | impbcor(imagename='3c391_ctm_spw0_multiscale.image',pbimage='3c391_ctm_spw0_multiscale.pb', |
300 + | outfile='3c391_ctm_spw0_multiscale.pbcorimage') |
301 + | |
302 + | mystat = imstat(imagename='3c391_ctm_spw0_multiscale.pbcorimage') |
303 + | |
304 + | # TESTS: compare values in return dictionary |
305 + | |
306 + | # location of image peak |
307 + | result, expected = mystat['maxpos'][:2], np.array([288,256]) |
308 + | msg = f"Expected image peak at pixels {expected}, got {result}, tolerance exact" |
309 + | print(msg); report.append(msg) |
310 + | self.assertTrue( np.all(result==expected), msg=msg ) |
311 + | |
312 + | # value of image peak |
313 + | result, expected, epsilon = mystat['max'][0], 0.15553903579711914, 0.04 |
314 + | msg = f"Expected pbcor image peak of {expected}, got {result}, tolerance {100*epsilon}%" |
315 + | print(msg); report.append(msg) |
316 + | self.assertAlmostEqual(result, expected, delta=expected*epsilon, msg=msg ) |
317 + | |
318 + | |
319 + | mystat = imstat(imagename='3c391_ctm_spw0_multiscale.residual', region=self.maskfile) |
320 + | |
321 + | # TESTS: compare values in return dictionary |
322 + | |
323 + | # max residual inside mask |
324 + | result, expected, epsilon = mystat['max'][0], 0.0010077510960400105, 0.10 |
325 + | msg = f"Expected peak residual of {expected}, got {result}, tolerance {100*epsilon}%" |
326 + | print(msg); report.append(msg) |
327 + | self.assertAlmostEqual(result, expected, delta=expected*epsilon, msg=msg ) |
328 + | |
329 + | # peak residual inside mask |
330 + | result, expected, epsilon = mystat['rms'][0], 0.000542584067811474, 0.10 |
331 + | msg = f"Expected rms residual of {expected}, got {result}, tolerance {100*epsilon}%" |
332 + | print(msg); report.append(msg) |
333 + | self.assertAlmostEqual(result, expected, delta=expected*epsilon, msg=msg ) |
334 + | |
335 + | test_dict[test_name]['report'] = '\n'.join(report) |
336 + | test_dict[test_name]['images'] = [] |
337 + | |
338 + | # main |
339 + | if __name__ == '__main__': |
340 + | unittest.main() |