source: branches/alma/python/asapmath.py@ 1719

Last change on this file since 1719 was 1693, checked in by Takeshi Nakazato, 15 years ago

New Development: No

JIRA Issue: Yes CAS-1908

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes

Module(s): Module Names change impacts.

Description: Describe your changes here...

Changed a tagging as the source type is tagged by using SRCTYPE, not
an extra string in SRCNAME. To do this, I have defined a selection method
by SRCTYPE in STSelector class. I have newly added python_SrcType.cpp
that defines a Python interface of SrcType enums which is defined
in atnf/PKSIO/SrcType.h.

Since I have added new file in the src directory, I have modified src/Makefile
to compile new file.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 39.6 KB
Line 
1from asap.scantable import scantable
2from asap import rcParams
3from asap import print_log
4from asap import selector
5from asap import asaplog
6from asap import asaplotgui
7
8def average_time(*args, **kwargs):
9 """
10 Return the (time) average of a scan or list of scans. [in channels only]
11 The cursor of the output scan is set to 0
12 Parameters:
13 one scan or comma separated scans or a list of scans
14 mask: an optional mask (only used for 'var' and 'tsys' weighting)
15 scanav: True averages each scan separately.
16 False (default) averages all scans together,
17 weight: Weighting scheme.
18 'none' (mean no weight)
19 'var' (1/var(spec) weighted)
20 'tsys' (1/Tsys**2 weighted)
21 'tint' (integration time weighted)
22 'tintsys' (Tint/Tsys**2)
23 'median' ( median averaging)
24 align: align the spectra in velocity before averaging. It takes
25 the time of the first spectrum in the first scantable
26 as reference time.
27 Example:
28 # return a time averaged scan from scana and scanb
29 # without using a mask
30 scanav = average_time(scana,scanb)
31 # or equivalent
32 # scanav = average_time([scana, scanb])
33 # return the (time) averaged scan, i.e. the average of
34 # all correlator cycles
35 scanav = average_time(scan, scanav=True)
36 """
37 scanav = False
38 if kwargs.has_key('scanav'):
39 scanav = kwargs.get('scanav')
40 weight = 'tint'
41 if kwargs.has_key('weight'):
42 weight = kwargs.get('weight')
43 mask = ()
44 if kwargs.has_key('mask'):
45 mask = kwargs.get('mask')
46 align = False
47 if kwargs.has_key('align'):
48 align = kwargs.get('align')
49 compel = False
50 if kwargs.has_key('compel'):
51 compel = kwargs.get('compel')
52 varlist = vars()
53 if isinstance(args[0],list):
54 lst = args[0]
55 elif isinstance(args[0],tuple):
56 lst = list(args[0])
57 else:
58 lst = list(args)
59
60 del varlist["kwargs"]
61 varlist["args"] = "%d scantables" % len(lst)
62 # need special formatting here for history...
63
64 from asap._asap import stmath
65 stm = stmath()
66 for s in lst:
67 if not isinstance(s,scantable):
68 msg = "Please give a list of scantables"
69 if rcParams['verbose']:
70 #print msg
71 asaplog.push(msg)
72 print_log('ERROR')
73 return
74 else:
75 raise TypeError(msg)
76 if scanav: scanav = "SCAN"
77 else: scanav = "NONE"
78 alignedlst = []
79 if align:
80 refepoch = lst[0].get_time(0)
81 for scan in lst:
82 alignedlst.append(scan.freq_align(refepoch,insitu=False))
83 else:
84 alignedlst = lst
85 if weight.upper() == 'MEDIAN':
86 # median doesn't support list of scantables - merge first
87 merged = None
88 if len(alignedlst) > 1:
89 merged = merge(alignedlst)
90 else:
91 merged = alignedlst[0]
92 s = scantable(stm._averagechannel(merged, 'MEDIAN', scanav))
93 del merged
94 else:
95 #s = scantable(stm._average(alignedlst, mask, weight.upper(), scanav))
96 s = scantable(stm._new_average(alignedlst, compel, mask, weight.upper(), scanav))
97 s._add_history("average_time",varlist)
98 print_log()
99 return s
100
101def quotient(source, reference, preserve=True):
102 """
103 Return the quotient of a 'source' (signal) scan and a 'reference' scan.
104 The reference can have just one scan, even if the signal has many. Otherwise
105 they must have the same number of scans.
106 The cursor of the output scan is set to 0
107 Parameters:
108 source: the 'on' scan
109 reference: the 'off' scan
110 preserve: you can preserve (default) the continuum or
111 remove it. The equations used are
112 preserve: Output = Toff * (on/off) - Toff
113 remove: Output = Toff * (on/off) - Ton
114 """
115 varlist = vars()
116 from asap._asap import stmath
117 stm = stmath()
118 stm._setinsitu(False)
119 s = scantable(stm._quotient(source, reference, preserve))
120 s._add_history("quotient",varlist)
121 print_log()
122 return s
123
124def dototalpower(calon, caloff, tcalval=0.0):
125 """
126 Do calibration for CAL on,off signals.
127 Adopted from GBTIDL dototalpower
128 Parameters:
129 calon: the 'cal on' subintegration
130 caloff: the 'cal off' subintegration
131 tcalval: user supplied Tcal value
132 """
133 varlist = vars()
134 from asap._asap import stmath
135 stm = stmath()
136 stm._setinsitu(False)
137 s = scantable(stm._dototalpower(calon, caloff, tcalval))
138 s._add_history("dototalpower",varlist)
139 print_log()
140 return s
141
142def dosigref(sig, ref, smooth, tsysval=0.0, tauval=0.0):
143 """
144 Calculate a quotient (sig-ref/ref * Tsys)
145 Adopted from GBTIDL dosigref
146 Parameters:
147 sig: on source data
148 ref: reference data
149 smooth: width of box car smoothing for reference
150 tsysval: user specified Tsys (scalar only)
151 tauval: user specified Tau (required if tsysval is set)
152 """
153 varlist = vars()
154 from asap._asap import stmath
155 stm = stmath()
156 stm._setinsitu(False)
157 s = scantable(stm._dosigref(sig, ref, smooth, tsysval, tauval))
158 s._add_history("dosigref",varlist)
159 print_log()
160 return s
161
162def calps(scantab, scannos, smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
163 """
164 Calibrate GBT position switched data
165 Adopted from GBTIDL getps
166 Currently calps identify the scans as position switched data if they
167 contain '_ps' in the source name. The data must contains 'CAL' signal
168 on/off in each integration. To identify 'CAL' on state, the word, 'calon'
169 need to be present in the source name field.
170 (GBT MS data reading process to scantable automatically append these
171 id names to the source names)
172
173 Parameters:
174 scantab: scantable
175 scannos: list of scan numbers
176 smooth: optional box smoothing order for the reference
177 (default is 1 = no smoothing)
178 tsysval: optional user specified Tsys (default is 0.0,
179 use Tsys in the data)
180 tauval: optional user specified Tau
181 tcalval: optional user specified Tcal (default is 0.0,
182 use Tcal value in the data)
183 """
184 varlist = vars()
185 # check for the appropriate data
186## s = scantab.get_scan('*_ps*')
187## if s is None:
188## msg = "The input data appear to contain no position-switch mode data."
189## if rcParams['verbose']:
190## #print msg
191## asaplog.push(msg)
192## print_log('ERROR')
193## return
194## else:
195## raise TypeError(msg)
196 s = scantab.copy()
197 from asap._asap import srctype
198 sel = selector()
199 sel.set_types( srctype.pson )
200 try:
201 scantab.set_selection( sel )
202 except Exception, e:
203 msg = "The input data appear to contain no position-switch mode data."
204 if rcParams['verbose']:
205 #print msg
206 asaplog.push(msg)
207 print_log('ERROR')
208 return
209 else:
210 raise TypeError(msg)
211 s.set_selection()
212 sel.reset()
213 ssub = s.get_scan(scannos)
214 if ssub is None:
215 msg = "No data was found with given scan numbers!"
216 if rcParams['verbose']:
217 #print msg
218 asaplog.push(msg)
219 print_log('ERROR')
220 return
221 else:
222 raise TypeError(msg)
223 #ssubon = ssub.get_scan('*calon')
224 #ssuboff = ssub.get_scan('*[^calon]')
225 sel.set_types( [srctype.poncal,srctype.poffcal] )
226 ssub.set_selection( sel )
227 ssubon = ssub.copy()
228 ssub.set_selection()
229 sel.reset()
230 sel.set_types( [srctype.pson,srctype.psoff] )
231 ssub.set_selection( sel )
232 ssuboff = ssub.copy()
233 ssub.set_selection()
234 sel.reset()
235 if ssubon.nrow() != ssuboff.nrow():
236 msg = "mismatch in numbers of CAL on/off scans. Cannot calibrate. Check the scan numbers."
237 if rcParams['verbose']:
238 #print msg
239 asaplog.push(msg)
240 print_log('ERROR')
241 return
242 else:
243 raise TypeError(msg)
244 cals = dototalpower(ssubon, ssuboff, tcalval)
245 #sig = cals.get_scan('*ps')
246 #ref = cals.get_scan('*psr')
247 sel.set_types( srctype.pson )
248 cals.set_selection( sel )
249 sig = cals.copy()
250 cals.set_selection()
251 sel.reset()
252 sel.set_types( srctype.psoff )
253 cals.set_selection( sel )
254 ref = cals.copy()
255 cals.set_selection()
256 sel.reset()
257 if sig.nscan() != ref.nscan():
258 msg = "mismatch in numbers of on/off scans. Cannot calibrate. Check the scan numbers."
259 if rcParams['verbose']:
260 #print msg
261 asaplog.push(msg)
262 print_log('ERROR')
263 return
264 else:
265 raise TypeError(msg)
266
267 #for user supplied Tsys
268 if tsysval>0.0:
269 if tauval<=0.0:
270 msg = "Need to supply a valid tau to use the supplied Tsys"
271 if rcParams['verbose']:
272 #print msg
273 asaplog.push(msg)
274 print_log('ERROR')
275 return
276 else:
277 raise TypeError(msg)
278 else:
279 sig.recalc_azel()
280 ref.recalc_azel()
281 #msg = "Use of user specified Tsys is not fully implemented yet."
282 #if rcParams['verbose']:
283 # print msg
284 # return
285 #else:
286 # raise TypeError(msg)
287 # use get_elevation to get elevation and
288 # calculate a scaling factor using the formula
289 # -> tsys use to dosigref
290
291 #ress = dosigref(sig, ref, smooth, tsysval)
292 ress = dosigref(sig, ref, smooth, tsysval, tauval)
293 ###
294 if verify:
295 # get data
296 import numpy
297 precal={}
298 postcal=[]
299 keys=['ps','ps_calon','psr','psr_calon']
300 types=[srctype.pson,srctype.poncal,srctype.psoff,srctype.poffcal]
301 ifnos=list(ssub.getifnos())
302 polnos=list(ssub.getpolnos())
303 sel=selector()
304 for i in range(2):
305 #ss=ssuboff.get_scan('*'+keys[2*i])
306 ll=[]
307 for j in range(len(ifnos)):
308 for k in range(len(polnos)):
309 sel.set_ifs(ifnos[j])
310 sel.set_polarizations(polnos[k])
311 sel.set_types(types[2*i])
312 try:
313 #ss.set_selection(sel)
314 ssuboff.set_selection(sel)
315 except:
316 continue
317 #ll.append(numpy.array(ss._getspectrum(0)))
318 ll.append(numpy.array(ssuboff._getspectrum(0)))
319 sel.reset()
320 ssuboff.set_selection()
321 precal[keys[2*i]]=ll
322 #del ss
323 #ss=ssubon.get_scan('*'+keys[2*i+1])
324 ll=[]
325 for j in range(len(ifnos)):
326 for k in range(len(polnos)):
327 sel.set_ifs(ifnos[j])
328 sel.set_polarizations(polnos[k])
329 sel.set_types(types[2*i+1])
330 try:
331 #ss.set_selection(sel)
332 ssubon.set_selection(sel)
333 except:
334 continue
335 #ll.append(numpy.array(ss._getspectrum(0)))
336 ll.append(numpy.array(ssubon._getspectrum(0)))
337 sel.reset()
338 ssubon.set_selection()
339 precal[keys[2*i+1]]=ll
340 #del ss
341 for j in range(len(ifnos)):
342 for k in range(len(polnos)):
343 sel.set_ifs(ifnos[j])
344 sel.set_polarizations(polnos[k])
345 try:
346 ress.set_selection(sel)
347 except:
348 continue
349 postcal.append(numpy.array(ress._getspectrum(0)))
350 sel.reset()
351 ress.set_selection()
352 del sel
353 # plot
354 print_log()
355 asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.')
356 print_log('WARN')
357 p=asaplotgui.asaplotgui()
358 #nr=min(6,len(ifnos)*len(polnos))
359 nr=len(ifnos)*len(polnos)
360 titles=[]
361 btics=[]
362 if nr<4:
363 p.set_panels(rows=nr,cols=2,nplots=2*nr,ganged=False)
364 for i in range(2*nr):
365 b=False
366 if i >= 2*nr-2:
367 b=True
368 btics.append(b)
369 elif nr==4:
370 p.set_panels(rows=2,cols=4,nplots=8,ganged=False)
371 for i in range(2*nr):
372 b=False
373 if i >= 2*nr-4:
374 b=True
375 btics.append(b)
376 elif nr<7:
377 p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
378 for i in range(2*nr):
379 if i >= 2*nr-4:
380 b=True
381 btics.append(b)
382 else:
383 print_log()
384 asaplog.push('Only first 6 [if,pol] pairs are plotted.')
385 print_log('WARN')
386 nr=6
387 for i in range(2*nr):
388 b=False
389 if i >= 2*nr-4:
390 b=True
391 btics.append(b)
392 p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
393 for i in range(nr):
394 p.subplot(2*i)
395 p.color=0
396 title='raw data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
397 titles.append(title)
398 #p.set_axes('title',title,fontsize=40)
399 ymin=1.0e100
400 ymax=-1.0e100
401 nchan=s.nchan()
402 edge=int(nchan*0.01)
403 for j in range(4):
404 spmin=min(precal[keys[j]][i][edge:nchan-edge])
405 spmax=max(precal[keys[j]][i][edge:nchan-edge])
406 ymin=min(ymin,spmin)
407 ymax=max(ymax,spmax)
408 for j in range(4):
409 if i==0:
410 p.set_line(label=keys[j])
411 else:
412 p.legend()
413 p.plot(precal[keys[j]][i])
414 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
415 if not btics[2*i]:
416 p.axes.set_xticks([])
417 p.subplot(2*i+1)
418 p.color=0
419 title='cal data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
420 titles.append(title)
421 #p.set_axes('title',title)
422 p.legend()
423 ymin=postcal[i][edge:nchan-edge].min()
424 ymax=postcal[i][edge:nchan-edge].max()
425 p.plot(postcal[i])
426 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
427 if not btics[2*i+1]:
428 p.axes.set_xticks([])
429 for i in range(2*nr):
430 p.subplot(i)
431 p.set_axes('title',titles[i],fontsize='medium')
432 x=raw_input('Accept calibration ([y]/n): ' )
433 if x.upper() == 'N':
434 p.unmap()
435 del p
436 return scabtab
437 p.unmap()
438 del p
439 ###
440 ress._add_history("calps", varlist)
441 print_log()
442 return ress
443
444def calnod(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
445 """
446 Do full (but a pair of scans at time) processing of GBT Nod data
447 calibration.
448 Adopted from GBTIDL's getnod
449 Parameters:
450 scantab: scantable
451 scannos: a pair of scan numbers, or the first scan number of the pair
452 smooth: box car smoothing order
453 tsysval: optional user specified Tsys value
454 tauval: optional user specified tau value (not implemented yet)
455 tcalval: optional user specified Tcal value
456 """
457 varlist = vars()
458 from asap._asap import stmath
459 from asap._asap import srctype
460 stm = stmath()
461 stm._setinsitu(False)
462
463 # check for the appropriate data
464## s = scantab.get_scan('*_nod*')
465## if s is None:
466## msg = "The input data appear to contain no Nod observing mode data."
467## if rcParams['verbose']:
468## #print msg
469## asaplog.push(msg)
470## print_log('ERROR')
471## return
472## else:
473## raise TypeError(msg)
474 s = scantab.copy()
475 sel = selector()
476 sel.set_types( srctype.nod )
477 try:
478 s.set_selection( sel )
479 except Exception, e:
480 msg = "The input data appear to contain no Nod observing mode data."
481 if rcParams['verbose']:
482 #print msg
483 asaplog.push(msg)
484 print_log('ERROR')
485 return
486 else:
487 raise TypeError(msg)
488 sel.reset()
489 del sel
490 del s
491
492 # need check correspondance of each beam with sig-ref ...
493 # check for timestamps, scan numbers, subscan id (not available in
494 # ASAP data format...). Assume 1st scan of the pair (beam 0 - sig
495 # and beam 1 - ref...)
496 # First scan number of paired scans or list of pairs of
497 # scan numbers (has to have even number of pairs.)
498
499 #data splitting
500 scan1no = scan2no = 0
501
502 if len(scannos)==1:
503 scan1no = scannos[0]
504 scan2no = scannos[0]+1
505 pairScans = [scan1no, scan2no]
506 else:
507 #if len(scannos)>2:
508 # msg = "calnod can only process a pair of nod scans at time."
509 # if rcParams['verbose']:
510 # print msg
511 # return
512 # else:
513 # raise TypeError(msg)
514 #
515 #if len(scannos)==2:
516 # scan1no = scannos[0]
517 # scan2no = scannos[1]
518 pairScans = list(scannos)
519
520 if tsysval>0.0:
521 if tauval<=0.0:
522 msg = "Need to supply a valid tau to use the supplied Tsys"
523 if rcParams['verbose']:
524 #print msg
525 asaplog.push(msg)
526 print_log('ERROR')
527 return
528 else:
529 raise TypeError(msg)
530 else:
531 scantab.recalc_azel()
532 resspec = scantable(stm._donod(scantab, pairScans, smooth, tsysval,tauval,tcalval))
533 ###
534 if verify:
535 # get data
536 import numpy
537 precal={}
538 postcal=[]
539 keys=['','_calon']
540 types=[srctype.nod,srctype.nodcal]
541 ifnos=list(scantab.getifnos())
542 polnos=list(scantab.getpolnos())
543 sel=selector()
544 ss = scantab.copy()
545 for i in range(2):
546 #ss=scantab.get_scan('*'+keys[i])
547 ll=[]
548 ll2=[]
549 for j in range(len(ifnos)):
550 for k in range(len(polnos)):
551 sel.set_ifs(ifnos[j])
552 sel.set_polarizations(polnos[k])
553 sel.set_scans(pairScans[0])
554 sel.set_types(types[i])
555 try:
556 ss.set_selection(sel)
557 except:
558 continue
559 ll.append(numpy.array(ss._getspectrum(0)))
560 sel.reset()
561 ss.set_selection()
562 sel.set_ifs(ifnos[j])
563 sel.set_polarizations(polnos[k])
564 sel.set_scans(pairScans[1])
565 sel.set_types(types[i])
566 try:
567 ss.set_selection(sel)
568 except:
569 ll.pop()
570 continue
571 ll2.append(numpy.array(ss._getspectrum(0)))
572 sel.reset()
573 ss.set_selection()
574 key='%s%s' %(pairScans[0],keys[i])
575 precal[key]=ll
576 key='%s%s' %(pairScans[1],keys[i])
577 precal[key]=ll2
578 #del ss
579 keys=precal.keys()
580 for j in range(len(ifnos)):
581 for k in range(len(polnos)):
582 sel.set_ifs(ifnos[j])
583 sel.set_polarizations(polnos[k])
584 sel.set_scans(pairScans[0])
585 try:
586 resspec.set_selection(sel)
587 except:
588 continue
589 postcal.append(numpy.array(resspec._getspectrum(0)))
590 sel.reset()
591 resspec.set_selection()
592 del sel
593 # plot
594 print_log()
595 asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.')
596 print_log('WARN')
597 p=asaplotgui.asaplotgui()
598 #nr=min(6,len(ifnos)*len(polnos))
599 nr=len(ifnos)*len(polnos)
600 titles=[]
601 btics=[]
602 if nr<4:
603 p.set_panels(rows=nr,cols=2,nplots=2*nr,ganged=False)
604 for i in range(2*nr):
605 b=False
606 if i >= 2*nr-2:
607 b=True
608 btics.append(b)
609 elif nr==4:
610 p.set_panels(rows=2,cols=4,nplots=8,ganged=False)
611 for i in range(2*nr):
612 b=False
613 if i >= 2*nr-4:
614 b=True
615 btics.append(b)
616 elif nr<7:
617 p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
618 for i in range(2*nr):
619 if i >= 2*nr-4:
620 b=True
621 btics.append(b)
622 else:
623 print_log()
624 asaplog.push('Only first 6 [if,pol] pairs are plotted.')
625 print_log('WARN')
626 nr=6
627 for i in range(2*nr):
628 b=False
629 if i >= 2*nr-4:
630 b=True
631 btics.append(b)
632 p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
633 for i in range(nr):
634 p.subplot(2*i)
635 p.color=0
636 title='raw data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
637 titles.append(title)
638 #p.set_axes('title',title,fontsize=40)
639 ymin=1.0e100
640 ymax=-1.0e100
641 nchan=scantab.nchan()
642 edge=int(nchan*0.01)
643 for j in range(4):
644 spmin=min(precal[keys[j]][i][edge:nchan-edge])
645 spmax=max(precal[keys[j]][i][edge:nchan-edge])
646 ymin=min(ymin,spmin)
647 ymax=max(ymax,spmax)
648 for j in range(4):
649 if i==0:
650 p.set_line(label=keys[j])
651 else:
652 p.legend()
653 p.plot(precal[keys[j]][i])
654 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
655 if not btics[2*i]:
656 p.axes.set_xticks([])
657 p.subplot(2*i+1)
658 p.color=0
659 title='cal data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
660 titles.append(title)
661 #p.set_axes('title',title)
662 p.legend()
663 ymin=postcal[i][edge:nchan-edge].min()
664 ymax=postcal[i][edge:nchan-edge].max()
665 p.plot(postcal[i])
666 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
667 if not btics[2*i+1]:
668 p.axes.set_xticks([])
669 for i in range(2*nr):
670 p.subplot(i)
671 p.set_axes('title',titles[i],fontsize='medium')
672 x=raw_input('Accept calibration ([y]/n): ' )
673 if x.upper() == 'N':
674 p.unmap()
675 del p
676 return scabtab
677 p.unmap()
678 del p
679 ###
680 resspec._add_history("calnod",varlist)
681 print_log()
682 return resspec
683
684def calfs(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
685 """
686 Calibrate GBT frequency switched data.
687 Adopted from GBTIDL getfs.
688 Currently calfs identify the scans as frequency switched data if they
689 contain '_fs' in the source name. The data must contains 'CAL' signal
690 on/off in each integration. To identify 'CAL' on state, the word, 'calon'
691 need to be present in the source name field.
692 (GBT MS data reading via scantable automatically append these
693 id names to the source names)
694
695 Parameters:
696 scantab: scantable
697 scannos: list of scan numbers
698 smooth: optional box smoothing order for the reference
699 (default is 1 = no smoothing)
700 tsysval: optional user specified Tsys (default is 0.0,
701 use Tsys in the data)
702 tauval: optional user specified Tau
703 """
704 varlist = vars()
705 from asap._asap import stmath
706 from asap._asap import srctype
707 stm = stmath()
708 stm._setinsitu(False)
709
710# check = scantab.get_scan('*_fs*')
711# if check is None:
712# msg = "The input data appear to contain no Nod observing mode data."
713# if rcParams['verbose']:
714# print msg
715# return
716# else:
717# raise TypeError(msg)
718 s = scantab.get_scan(scannos)
719 del scantab
720
721 resspec = scantable(stm._dofs(s, scannos, smooth, tsysval,tauval,tcalval))
722 ###
723 if verify:
724 # get data
725 ssub = s.get_scan(scannos)
726 #ssubon = ssub.get_scan('*calon')
727 #ssuboff = ssub.get_scan('*[^calon]')
728 sel = selector()
729 sel.set_types( [srctype.foncal,srctype.foffcal] )
730 ssub.set_selection( sel )
731 ssubon = ssub.copy()
732 ssub.set_selection()
733 sel.reset()
734 sel.set_types( [srctype.fson,srctype.fsoff] )
735 ssub.set_selection( sel )
736 ssuboff = ssub.copy()
737 ssub.set_selection()
738 sel.reset()
739 import numpy
740 precal={}
741 postcal=[]
742 keys=['fs','fs_calon','fsr','fsr_calon']
743 types=[srctype.fson,srctype.foncal,srctype.fsoff,srctype.foffcal]
744 ifnos=list(ssub.getifnos())
745 polnos=list(ssub.getpolnos())
746 for i in range(2):
747 #ss=ssuboff.get_scan('*'+keys[2*i])
748 ll=[]
749 for j in range(len(ifnos)):
750 for k in range(len(polnos)):
751 sel.set_ifs(ifnos[j])
752 sel.set_polarizations(polnos[k])
753 sel.set_types(types[2*i])
754 try:
755 #ss.set_selection(sel)
756 ssuboff.set_selection(sel)
757 except:
758 continue
759 ll.append(numpy.array(ss._getspectrum(0)))
760 sel.reset()
761 #ss.set_selection()
762 ssuboff.set_selection()
763 precal[keys[2*i]]=ll
764 #del ss
765 #ss=ssubon.get_scan('*'+keys[2*i+1])
766 ll=[]
767 for j in range(len(ifnos)):
768 for k in range(len(polnos)):
769 sel.set_ifs(ifnos[j])
770 sel.set_polarizations(polnos[k])
771 sel.set_types(types[2*i+1])
772 try:
773 #ss.set_selection(sel)
774 ssubon.set_selection(sel)
775 except:
776 continue
777 ll.append(numpy.array(ss._getspectrum(0)))
778 sel.reset()
779 #ss.set_selection()
780 ssubon.set_selection()
781 precal[keys[2*i+1]]=ll
782 #del ss
783 #sig=resspec.get_scan('*_fs')
784 #ref=resspec.get_scan('*_fsr')
785 sel.set_types( srctype.fson )
786 resspec.set_selection( sel )
787 sig=resspec.copy()
788 resspec.set_selection()
789 sel.reset()
790 sel.set_type( srctype.fsoff )
791 resspec.set_selection( sel )
792 ref=resspec.copy()
793 resspec.set_selection()
794 sel.reset()
795 for k in range(len(polnos)):
796 for j in range(len(ifnos)):
797 sel.set_ifs(ifnos[j])
798 sel.set_polarizations(polnos[k])
799 try:
800 sig.set_selection(sel)
801 postcal.append(numpy.array(sig._getspectrum(0)))
802 except:
803 ref.set_selection(sel)
804 postcal.append(numpy.array(ref._getspectrum(0)))
805 sel.reset()
806 resspec.set_selection()
807 del sel
808 # plot
809 print_log()
810 asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.')
811 print_log('WARN')
812 p=asaplotgui.asaplotgui()
813 #nr=min(6,len(ifnos)*len(polnos))
814 nr=len(ifnos)/2*len(polnos)
815 titles=[]
816 btics=[]
817 if nr>3:
818 print_log()
819 asaplog.push('Only first 3 [if,pol] pairs are plotted.')
820 print_log('WARN')
821 nr=3
822 p.set_panels(rows=nr,cols=3,nplots=3*nr,ganged=False)
823 for i in range(3*nr):
824 b=False
825 if i >= 3*nr-3:
826 b=True
827 btics.append(b)
828 for i in range(nr):
829 p.subplot(3*i)
830 p.color=0
831 title='raw data IF%s,%s POL%s' % (ifnos[2*int(i/len(polnos))],ifnos[2*int(i/len(polnos))+1],polnos[i%len(polnos)])
832 titles.append(title)
833 #p.set_axes('title',title,fontsize=40)
834 ymin=1.0e100
835 ymax=-1.0e100
836 nchan=s.nchan()
837 edge=int(nchan*0.01)
838 for j in range(4):
839 spmin=min(precal[keys[j]][i][edge:nchan-edge])
840 spmax=max(precal[keys[j]][i][edge:nchan-edge])
841 ymin=min(ymin,spmin)
842 ymax=max(ymax,spmax)
843 for j in range(4):
844 if i==0:
845 p.set_line(label=keys[j])
846 else:
847 p.legend()
848 p.plot(precal[keys[j]][i])
849 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
850 if not btics[3*i]:
851 p.axes.set_xticks([])
852 p.subplot(3*i+1)
853 p.color=0
854 title='sig data IF%s POL%s' % (ifnos[2*int(i/len(polnos))],polnos[i%len(polnos)])
855 titles.append(title)
856 #p.set_axes('title',title)
857 p.legend()
858 ymin=postcal[2*i][edge:nchan-edge].min()
859 ymax=postcal[2*i][edge:nchan-edge].max()
860 p.plot(postcal[2*i])
861 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
862 if not btics[3*i+1]:
863 p.axes.set_xticks([])
864 p.subplot(3*i+2)
865 p.color=0
866 title='ref data IF%s POL%s' % (ifnos[2*int(i/len(polnos))+1],polnos[i%len(polnos)])
867 titles.append(title)
868 #p.set_axes('title',title)
869 p.legend()
870 ymin=postcal[2*i+1][edge:nchan-edge].min()
871 ymax=postcal[2*i+1][edge:nchan-edge].max()
872 p.plot(postcal[2*i+1])
873 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
874 if not btics[3*i+2]:
875 p.axes.set_xticks([])
876 for i in range(3*nr):
877 p.subplot(i)
878 p.set_axes('title',titles[i],fontsize='medium')
879 x=raw_input('Accept calibration ([y]/n): ' )
880 if x.upper() == 'N':
881 p.unmap()
882 del p
883 return scabtab
884 p.unmap()
885 del p
886 ###
887 resspec._add_history("calfs",varlist)
888 print_log()
889 return resspec
890
891def simple_math(left, right, op='add', tsys=True):
892 """
893 Apply simple mathematical binary operations to two
894 scan tables, returning the result in a new scan table.
895 The operation is applied to both the correlations and the TSys data
896 The cursor of the output scan is set to 0
897 Parameters:
898 left: the 'left' scan
899 right: the 'right' scan
900 op: the operation: 'add' (default), 'sub', 'mul', 'div'
901 tsys: if True (default) then apply the operation to Tsys
902 as well as the data
903 """
904 #print "simple_math is deprecated use +=/* instead."
905 asaplog.push( "simple_math is deprecated use +=/* instead." )
906 print_log('WARN')
907
908def merge(*args):
909 """
910 Merge a list of scanatables, or comma-sperated scantables into one
911 scnatble.
912 Parameters:
913 A list [scan1, scan2] or scan1, scan2.
914 Example:
915 myscans = [scan1, scan2]
916 allscans = merge(myscans)
917 # or equivalent
918 sameallscans = merge(scan1, scan2)
919 """
920 varlist = vars()
921 if isinstance(args[0],list):
922 lst = tuple(args[0])
923 elif isinstance(args[0],tuple):
924 lst = args[0]
925 else:
926 lst = tuple(args)
927 varlist["args"] = "%d scantables" % len(lst)
928 # need special formatting her for history...
929 from asap._asap import stmath
930 stm = stmath()
931 for s in lst:
932 if not isinstance(s,scantable):
933 msg = "Please give a list of scantables"
934 if rcParams['verbose']:
935 #print msg
936 asaplog.push(msg)
937 print_log('ERROR')
938 return
939 else:
940 raise TypeError(msg)
941 s = scantable(stm._merge(lst))
942 s._add_history("merge", varlist)
943 print_log()
944 return s
945
946def calibrate( scantab, scannos=[], calmode='none', verify=None ):
947 """
948 Calibrate data.
949
950 Parameters:
951 scantab: scantable
952 scannos: list of scan number
953 calmode: calibration mode
954 verify: verify calibration
955 """
956 antname = scantab.get_antennaname()
957 if ( calmode == 'nod' ):
958 asaplog.push( 'Calibrating nod data.' )
959 print_log()
960 scal = calnod( scantab, scannos=scannos, verify=verify )
961 elif ( calmode == 'quotient' ):
962 asaplog.push( 'Calibrating using quotient.' )
963 print_log()
964 scal = scantab.auto_quotient( verify=verify )
965 elif ( calmode == 'ps' ):
966 asaplog.push( 'Calibrating %s position-switched data.' % antname )
967 print_log()
968 if ( antname.find( 'APEX' ) != -1 ):
969 scal = apexcal( scantab, scannos, calmode, verify )
970 elif ( antname.find( 'ALMA' ) != -1 or antname.find( 'OSF' ) != -1 ):
971 scal = almacal( scantab, scannos, calmode, verify )
972 else:
973 scal = calps( scantab, scannos=scannos, verify=verify )
974 elif ( calmode == 'fs' or calmode == 'fsotf' ):
975 asaplog.push( 'Calibrating %s frequency-switched data.' % antname )
976 print_log()
977 if ( antname.find( 'APEX' ) != -1 ):
978 scal = apexcal( scantab, scannos, calmode, verify )
979 elif ( antname.find( 'ALMA' ) != -1 or antname.find( 'OSF' ) != -1 ):
980 scal = almacal( scantab, scannos, calmode, verify )
981 else:
982 scal = calfs( scantab, scannos=scannos, verify=verify )
983 elif ( calmode == 'otf' ):
984 asaplog.push( 'Calibrating %s On-The-Fly data.' % antname )
985 print_log()
986 scal = almacal( scantab, scannos, calmode, verify )
987 else:
988 asaplog.push( 'No calibration.' )
989 scal = scantab.copy()
990
991 return scal
992
993def apexcal( scantab, scannos=[], calmode='none', verify=False ):
994 """
995 Calibrate APEX data
996
997 Parameters:
998 scantab: scantable
999 scannos: list of scan number
1000 calmode: calibration mode
1001
1002 verify: verify calibration
1003 """
1004 from asap._asap import stmath
1005 stm = stmath()
1006 antname = scantab.get_antennaname()
1007 ssub = scantab.get_scan( scannos )
1008 scal = scantable( stm.cwcal( ssub, calmode, antname ) )
1009 return scal
1010
1011def almacal( scantab, scannos=[], calmode='none', verify=False ):
1012 """
1013 Calibrate ALMA data
1014
1015 Parameters:
1016 scantab: scantable
1017 scannos: list of scan number
1018 calmode: calibration mode
1019
1020 verify: verify calibration
1021 """
1022 from asap._asap import stmath
1023 stm = stmath()
1024 ssub = scantab.get_scan( scannos )
1025 scal = scantable( stm.almacal( ssub, calmode ) )
1026 return scal
1027
1028def splitant(filename, outprefix='',overwrite=False):
1029 """
1030 Split Measurement set by antenna name, save data as a scantables,
1031 and return a list of filename.
1032 Notice this method can only be available from CASA.
1033 Prameter
1034 filename: the name of Measurement set to be read.
1035 outprefix: the prefix of output scantable name.
1036 the names of output scantable will be
1037 outprefix.antenna1, outprefix.antenna2, ....
1038 If not specified, outprefix = filename is assumed.
1039 overwrite If the file should be overwritten if it exists.
1040 The default False is to return with warning
1041 without writing the output. USE WITH CARE.
1042
1043 """
1044 # Import the table toolkit from CASA
1045 try:
1046 import casac
1047 except ImportError:
1048 if rcParams['verbose']:
1049 #print "failed to load casa"
1050 print_log()
1051 asaplog.push("failed to load casa")
1052 print_log('ERROR')
1053 else: raise
1054 return False
1055 try:
1056 tbtool = casac.homefinder.find_home_by_name('tableHome')
1057 tb = tbtool.create()
1058 tb2 = tbtool.create()
1059 except:
1060 if rcParams['verbose']:
1061 #print "failed to load a table tool:\n", e
1062 print_log()
1063 asaplog.push("failed to load table tool")
1064 print_log('ERROR')
1065 else: raise
1066 return False
1067 # Check the input filename
1068 if isinstance(filename, str):
1069 import os.path
1070 filename = os.path.expandvars(filename)
1071 filename = os.path.expanduser(filename)
1072 if not os.path.exists(filename):
1073 s = "File '%s' not found." % (filename)
1074 if rcParams['verbose']:
1075 print_log()
1076 asaplog.push(s)
1077 print_log('ERROR')
1078 return
1079 raise IOError(s)
1080 # check if input file is MS
1081 if not os.path.isdir(filename) \
1082 or not os.path.exists(filename+'/ANTENNA') \
1083 or not os.path.exists(filename+'/table.f1'):
1084 s = "File '%s' is not a Measurement set." % (filename)
1085 if rcParams['verbose']:
1086 print_log()
1087 asaplog.push(s)
1088 print_log('ERROR')
1089 return
1090 raise IOError(s)
1091 else:
1092 s = "The filename should be string. "
1093 if rcParams['verbose']:
1094 print_log()
1095 asaplog.push(s)
1096 print_log('ERROR')
1097 return
1098 raise TypeError(s)
1099 # Check out put file name
1100 outname=''
1101 if len(outprefix) > 0: prefix=outprefix+'.'
1102 else:
1103 prefix=filename.rstrip('/')
1104 # Now do the actual splitting.
1105 outfiles=[]
1106 tb.open(tablename=filename+'/ANTENNA',nomodify=True)
1107 nant=tb.nrows()
1108 antnames=tb.getcol('NAME',0,nant,1)
1109 antpos=tb.getcol('POSITION',0,nant,1).transpose()
1110 tb.close()
1111 tb.open(tablename=filename,nomodify=True)
1112 ant1=tb.getcol('ANTENNA1',0,-1,1)
1113 tb.close()
1114 for antid in set(ant1):
1115 scan=scantable(filename,average=False,getpt=True,antenna=int(antid))
1116 outname=prefix+antnames[antid]+'.asap'
1117 scan.save(outname,format='ASAP',overwrite=overwrite)
1118 del scan
1119 outfiles.append(outname)
1120 del tb, tb2
1121 return outfiles
1122
1123def _array2dOp( scan, value, mode="ADD", tsys=False ):
1124 """
1125 This function is workaround on the basic operation of scantable
1126 with 2 dimensional float list.
1127
1128 scan: scantable operand
1129 value: float list operand
1130 mode: operation mode (ADD, SUB, MUL, DIV)
1131 tsys: if True, operate tsys as well
1132 """
1133 nrow = scan.nrow()
1134 s = None
1135 if len( value ) == 1:
1136 from asap._asap import stmath
1137 stm = stmath()
1138 s = scantable( stm._arrayop( scan.copy(), value[0], mode, tsys ) )
1139 del stm
1140 elif len( value ) != nrow:
1141 asaplog.push( 'len(value) must be 1 or conform to scan.nrow()' )
1142 print_log( 'ERROR' )
1143 else:
1144 from asap._asap import stmath
1145 stm = stmath()
1146 # insitu must be True
1147 stm._setinsitu( True )
1148 s = scan.copy()
1149 sel = selector()
1150 for irow in range( nrow ):
1151 sel.set_rows( irow )
1152 s.set_selection( sel )
1153 if len( value[irow] ) == 1:
1154 stm._unaryop( s, value[irow][0], mode, tsys )
1155 else:
1156 stm._arrayop( s, value[irow], mode, tsys, 'channel' )
1157 s.set_selection()
1158 sel.reset()
1159 del sel
1160 del stm
1161 return s
1162
1163
1164
Note: See TracBrowser for help on using the repository browser.