source: trunk/python/asapmath.py@ 1837

Last change on this file since 1837 was 1827, checked in by Malte Marquarding, 14 years ago

fixed type in import name

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