source: branches/hpc34/python/asapmath.py@ 2849

Last change on this file since 2849 was 2580, checked in by ShinnosukeKawakami, 12 years ago

hpc33 merged asap-trunk

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