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

Last change on this file since 1892 was 1782, checked in by Kana Sugimoto, 14 years ago

New Development: Yes

JIRA Issue: No (merge)

Ready for Test: Yes

Interface Changes: No

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s): single dish package

Description:

Merged some unimplemented developments from the trunk.
[Details]

  1. Restored PKSmsg.cc, PKSmsg.h to alma branch from trunk
  2. Updated pks_matchs.cc and src/SConscript to ones in trunk
  3. uncommented "@print_log_dec"
  4. implemented bug fixes to add(), sub(), mul(), and div() in scantable.py
  5. implemented get_region() and get_point() in asaplotbase.py
  6. activated interactive keyword for text(), arrow(), annotate(), axvline(), axhline(), axvspan(), and axhspan() in asapplotter.py


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