source: trunk/python/asapmath.py@ 2283

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

New Development: No

JIRA Issue: No (a fix)

Ready for Test: Yes

Interface Changes: added a new method new_asaplot() to module asapplotter which returns asaplot instance based on your backend selection.

What Interface Changed:

Test Programs: sdaverage, sdsmooth with verify=True and plotlevel > 0

Put in Release Notes: No

Module(s):

Description: proper handling of plotting in non-TkAgg backend.


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