source: branches/plotter2/python/asapmath.py@ 3076

Last change on this file since 3076 was 2818, checked in by Malte Marquarding, 11 years ago

Issue #291: added scantable.set_sourcename to overwrite sourcename to allow freq_align to work. Exposed 'SOURCE' averaging to python api. Added some logging to freq_align refering to the sources it aligns to

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 37.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
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,
93 weight.upper(), scanav))
94 s._add_history("average_time",varlist)
95
96 return s
97
98@asaplog_post_dec
99def quotient(source, reference, preserve=True):
100 """
101 Return the quotient of a 'source' (signal) scan and a 'reference' scan.
102 The reference can have just one scan, even if the signal has many. Otherwise
103 they must have the same number of scans.
104 The cursor of the output scan is set to 0
105 Parameters:
106 source: the 'on' scan
107 reference: the 'off' scan
108 preserve: you can preserve (default) the continuum or
109 remove it. The equations used are
110 preserve: Output = Toff * (on/off) - Toff
111 remove: Output = Toff * (on/off) - Ton
112 """
113 varlist = vars()
114 from asap._asap import stmath
115 stm = stmath()
116 stm._setinsitu(False)
117 s = scantable(stm._quotient(source, reference, preserve))
118 s._add_history("quotient",varlist)
119 return s
120
121@asaplog_post_dec
122def dototalpower(calon, caloff, tcalval=0.0):
123 """
124 Do calibration for CAL on,off signals.
125 Adopted from GBTIDL dototalpower
126 Parameters:
127 calon: the 'cal on' subintegration
128 caloff: the 'cal off' subintegration
129 tcalval: user supplied Tcal value
130 """
131 varlist = vars()
132 from asap._asap import stmath
133 stm = stmath()
134 stm._setinsitu(False)
135 s = scantable(stm._dototalpower(calon, caloff, tcalval))
136 s._add_history("dototalpower",varlist)
137 return s
138
139@asaplog_post_dec
140def dosigref(sig, ref, smooth, tsysval=0.0, tauval=0.0):
141 """
142 Calculate a quotient (sig-ref/ref * Tsys)
143 Adopted from GBTIDL dosigref
144 Parameters:
145 sig: on source data
146 ref: reference data
147 smooth: width of box car smoothing for reference
148 tsysval: user specified Tsys (scalar only)
149 tauval: user specified Tau (required if tsysval is set)
150 """
151 varlist = vars()
152 from asap._asap import stmath
153 stm = stmath()
154 stm._setinsitu(False)
155 s = scantable(stm._dosigref(sig, ref, smooth, tsysval, tauval))
156 s._add_history("dosigref",varlist)
157 return s
158
159@asaplog_post_dec
160def calps(scantab, scannos, smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
161 """
162 Calibrate GBT position switched data
163 Adopted from GBTIDL getps
164 Currently calps identify the scans as position switched data if source
165 type enum is pson or psoff. The data must contains 'CAL' signal
166 on/off in each integration. To identify 'CAL' on state, the source type
167 enum of poncal and poffcal need to be present.
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=new_asaplot()
315 rcp('lines', linewidth=1)
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(ifnos[int(i/len(polnos))])
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=new_asaplot()
535 rcp('lines', linewidth=1)
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(ifnos[int(i/len(polnos))])
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.
630
631 Parameters:
632 scantab: scantable
633 scannos: list of scan numbers
634 smooth: optional box smoothing order for the reference
635 (default is 1 = no smoothing)
636 tsysval: optional user specified Tsys (default is 0.0,
637 use Tsys in the data)
638 tauval: optional user specified Tau
639 verify: Verify calibration if true
640 """
641 varlist = vars()
642 from asap._asap import stmath
643 from asap._asap import srctype
644 stm = stmath()
645 stm._setinsitu(False)
646
647# check = scantab.get_scan('*_fs*')
648# if check is None:
649# msg = "The input data appear to contain no Nod observing mode data."
650# raise TypeError(msg)
651 s = scantab.get_scan(scannos)
652 del scantab
653
654 resspec = scantable(stm._dofs(s, scannos, smooth, tsysval,tauval,tcalval))
655 ###
656 if verify:
657 # get data
658 ssub = s.get_scan(scannos)
659 #ssubon = ssub.get_scan('*calon')
660 #ssuboff = ssub.get_scan('*[^calon]')
661 sel = selector()
662 sel.set_types( [srctype.foncal,srctype.foffcal] )
663 ssub.set_selection( sel )
664 ssubon = ssub.copy()
665 ssub.set_selection()
666 sel.reset()
667 sel.set_types( [srctype.fson,srctype.fsoff] )
668 ssub.set_selection( sel )
669 ssuboff = ssub.copy()
670 ssub.set_selection()
671 sel.reset()
672 import numpy
673 precal={}
674 postcal=[]
675 keys=['fs','fs_calon','fsr','fsr_calon']
676 types=[srctype.fson,srctype.foncal,srctype.fsoff,srctype.foffcal]
677 ifnos=list(ssub.getifnos())
678 polnos=list(ssub.getpolnos())
679 for i in range(2):
680 #ss=ssuboff.get_scan('*'+keys[2*i])
681 ll=[]
682 for j in range(len(ifnos)):
683 for k in range(len(polnos)):
684 sel.set_ifs(ifnos[j])
685 sel.set_polarizations(polnos[k])
686 sel.set_types(types[2*i])
687 try:
688 #ss.set_selection(sel)
689 ssuboff.set_selection(sel)
690 except:
691 continue
692 ll.append(numpy.array(ss._getspectrum(0)))
693 sel.reset()
694 #ss.set_selection()
695 ssuboff.set_selection()
696 precal[keys[2*i]]=ll
697 #del ss
698 #ss=ssubon.get_scan('*'+keys[2*i+1])
699 ll=[]
700 for j in range(len(ifnos)):
701 for k in range(len(polnos)):
702 sel.set_ifs(ifnos[j])
703 sel.set_polarizations(polnos[k])
704 sel.set_types(types[2*i+1])
705 try:
706 #ss.set_selection(sel)
707 ssubon.set_selection(sel)
708 except:
709 continue
710 ll.append(numpy.array(ss._getspectrum(0)))
711 sel.reset()
712 #ss.set_selection()
713 ssubon.set_selection()
714 precal[keys[2*i+1]]=ll
715 #del ss
716 #sig=resspec.get_scan('*_fs')
717 #ref=resspec.get_scan('*_fsr')
718 sel.set_types( srctype.fson )
719 resspec.set_selection( sel )
720 sig=resspec.copy()
721 resspec.set_selection()
722 sel.reset()
723 sel.set_type( srctype.fsoff )
724 resspec.set_selection( sel )
725 ref=resspec.copy()
726 resspec.set_selection()
727 sel.reset()
728 for k in range(len(polnos)):
729 for j in range(len(ifnos)):
730 sel.set_ifs(ifnos[j])
731 sel.set_polarizations(polnos[k])
732 try:
733 sig.set_selection(sel)
734 postcal.append(numpy.array(sig._getspectrum(0)))
735 except:
736 ref.set_selection(sel)
737 postcal.append(numpy.array(ref._getspectrum(0)))
738 sel.reset()
739 resspec.set_selection()
740 del sel
741 # plot
742 asaplog.post()
743 asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.')
744 asaplog.post('WARN')
745 p=new_asaplot()
746 rcp('lines', linewidth=1)
747 #nr=min(6,len(ifnos)*len(polnos))
748 nr=len(ifnos)/2*len(polnos)
749 titles=[]
750 btics=[]
751 if nr>3:
752 asaplog.post()
753 asaplog.push('Only first 3 [if,pol] pairs are plotted.')
754 asaplog.post('WARN')
755 nr=3
756 p.set_panels(rows=nr,cols=3,nplots=3*nr,ganged=False)
757 for i in range(3*nr):
758 b=False
759 if i >= 3*nr-3:
760 b=True
761 btics.append(b)
762 for i in range(nr):
763 p.subplot(3*i)
764 p.color=0
765 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)])
766 titles.append(title)
767 #p.set_axes('title',title,fontsize=40)
768 ymin=1.0e100
769 ymax=-1.0e100
770 nchan=s.nchan(ifnos[2*int(i/len(polnos))])
771 edge=int(nchan*0.01)
772 for j in range(4):
773 spmin=min(precal[keys[j]][i][edge:nchan-edge])
774 spmax=max(precal[keys[j]][i][edge:nchan-edge])
775 ymin=min(ymin,spmin)
776 ymax=max(ymax,spmax)
777 for j in range(4):
778 if i==0:
779 p.set_line(label=keys[j])
780 else:
781 p.legend()
782 p.plot(precal[keys[j]][i])
783 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
784 if not btics[3*i]:
785 p.axes.set_xticks([])
786 p.subplot(3*i+1)
787 p.color=0
788 title='sig data IF%s POL%s' % (ifnos[2*int(i/len(polnos))],polnos[i%len(polnos)])
789 titles.append(title)
790 #p.set_axes('title',title)
791 p.legend()
792 ymin=postcal[2*i][edge:nchan-edge].min()
793 ymax=postcal[2*i][edge:nchan-edge].max()
794 p.plot(postcal[2*i])
795 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
796 if not btics[3*i+1]:
797 p.axes.set_xticks([])
798 p.subplot(3*i+2)
799 p.color=0
800 title='ref data IF%s POL%s' % (ifnos[2*int(i/len(polnos))+1],polnos[i%len(polnos)])
801 titles.append(title)
802 #p.set_axes('title',title)
803 p.legend()
804 ymin=postcal[2*i+1][edge:nchan-edge].min()
805 ymax=postcal[2*i+1][edge:nchan-edge].max()
806 p.plot(postcal[2*i+1])
807 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
808 if not btics[3*i+2]:
809 p.axes.set_xticks([])
810 for i in range(3*nr):
811 p.subplot(i)
812 p.set_axes('title',titles[i],fontsize='medium')
813 x=raw_input('Accept calibration ([y]/n): ' )
814 if x.upper() == 'N':
815 p.quit()
816 del p
817 return scabtab
818 p.quit()
819 del p
820 ###
821 resspec._add_history("calfs",varlist)
822 return resspec
823
824@asaplog_post_dec
825def merge(*args):
826 """
827 Merge a list of scanatables, or comma-sperated scantables into one
828 scnatble.
829 Parameters:
830 A list [scan1, scan2] or scan1, scan2.
831 Example:
832 myscans = [scan1, scan2]
833 allscans = merge(myscans)
834 # or equivalent
835 sameallscans = merge(scan1, scan2)
836 """
837 varlist = vars()
838 if isinstance(args[0],list):
839 lst = tuple(args[0])
840 elif isinstance(args[0],tuple):
841 lst = args[0]
842 else:
843 lst = tuple(args)
844 varlist["args"] = "%d scantables" % len(lst)
845 # need special formatting her for history...
846 from asap._asap import stmath
847 stm = stmath()
848 for s in lst:
849 if not isinstance(s,scantable):
850 msg = "Please give a list of scantables"
851 raise TypeError(msg)
852 s = scantable(stm._merge(lst))
853 s._add_history("merge", varlist)
854 return s
855
856@asaplog_post_dec
857def calibrate( scantab, scannos=[], calmode='none', verify=None ):
858 """
859 Calibrate data.
860
861 Parameters:
862 scantab: scantable
863 scannos: list of scan number
864 calmode: calibration mode
865 verify: verify calibration
866 """
867 import re
868 antname = scantab.get_antennaname()
869 if ( calmode == 'nod' ):
870 asaplog.push( 'Calibrating nod data.' )
871 scal = calnod( scantab, scannos=scannos, verify=verify )
872 elif ( calmode == 'quotient' ):
873 asaplog.push( 'Calibrating using quotient.' )
874 scal = scantab.auto_quotient( verify=verify )
875 elif ( calmode == 'ps' ):
876 asaplog.push( 'Calibrating %s position-switched data.' % antname )
877 if ( antname.find( 'APEX' ) != -1 ):
878 scal = apexcal( scantab, scannos, calmode, verify )
879 elif ( antname.find( 'ALMA' ) != -1 or antname.find( 'OSF' ) != -1
880 or re.match('DV[0-9][0-9]$',antname) is not None
881 or re.match('PM[0-9][0-9]$',antname) is not None
882 or re.match('CM[0-9][0-9]$',antname) is not None
883 or re.match('DA[0-9][0-9]$',antname) is not None ):
884 scal = almacal( scantab, scannos, calmode, verify )
885 else:
886 scal = calps( scantab, scannos=scannos, verify=verify )
887 elif ( calmode == 'fs' or calmode == 'fsotf' ):
888 asaplog.push( 'Calibrating %s frequency-switched data.' % antname )
889 if ( antname.find( 'APEX' ) != -1 ):
890 scal = apexcal( scantab, scannos, calmode, verify )
891 elif ( antname.find( 'ALMA' ) != -1 or antname.find( 'OSF' ) != -1 ):
892 scal = almacal( scantab, scannos, calmode, verify )
893 else:
894 scal = calfs( scantab, scannos=scannos, verify=verify )
895 elif ( calmode == 'otf' ):
896 asaplog.push( 'Calibrating %s On-The-Fly data.' % antname )
897 scal = almacal( scantab, scannos, calmode, verify )
898 else:
899 asaplog.push( 'No calibration.' )
900 scal = scantab.copy()
901
902 return scal
903
904def apexcal( scantab, scannos=[], calmode='none', verify=False ):
905 """
906 Calibrate APEX data
907
908 Parameters:
909 scantab: scantable
910 scannos: list of scan number
911 calmode: calibration mode
912
913 verify: verify calibration
914 """
915 from asap._asap import stmath
916 stm = stmath()
917 antname = scantab.get_antennaname()
918 selection=selector()
919 selection.set_scans(scannos)
920 orig = scantab.get_selection()
921 scantab.set_selection(orig+selection)
922## ssub = scantab.get_scan( scannos )
923## scal = scantable( stm.cwcal( ssub, calmode, antname ) )
924 scal = scantable( stm.cwcal( scantab, calmode, antname ) )
925 scantab.set_selection(orig)
926 return scal
927
928def almacal( scantab, scannos=[], calmode='none', verify=False ):
929 """
930 Calibrate ALMA data
931
932 Parameters:
933 scantab: scantable
934 scannos: list of scan number
935 calmode: calibration mode
936
937 verify: verify calibration
938 """
939 from asap._asap import stmath
940 stm = stmath()
941 selection=selector()
942 selection.set_scans(scannos)
943 orig = scantab.get_selection()
944 scantab.set_selection(orig+selection)
945## ssub = scantab.get_scan( scannos )
946## scal = scantable( stm.almacal( ssub, calmode ) )
947 scal = scantable( stm.almacal( scantab, calmode ) )
948 scantab.set_selection(orig)
949 return scal
950
951@asaplog_post_dec
952def splitant(filename, outprefix='',overwrite=False, getpt=True, freq_tolsr=False):
953 """
954 Split Measurement set by antenna name, save data as a scantables,
955 and return a list of filename.
956 Notice this method can only be available from CASA.
957 Prameter
958 filename: the name of Measurement set to be read.
959 outprefix: the prefix of output scantable name.
960 the names of output scantable will be
961 outprefix.antenna1, outprefix.antenna2, ....
962 If not specified, outprefix = filename is assumed.
963 overwrite If the file should be overwritten if it exists.
964 The default False is to return with warning
965 without writing the output. USE WITH CARE.
966 getpt Whether to import direction from MS/POINTING
967 table or not. Default is True (import direction).
968 freq_tolsr Whether to convert frequency frame information
969 to LSRK or not. Default is False (import as is).
970
971 """
972 # Import the table toolkit from CASA
973 from casac import casac
974 from asap.scantable import is_ms
975 tb = casac.table()
976 # Check the input filename
977 if isinstance(filename, str):
978 import os.path
979 filename = os.path.expandvars(filename)
980 filename = os.path.expanduser(filename)
981 if not os.path.exists(filename):
982 s = "File '%s' not found." % (filename)
983 raise IOError(s)
984 # check if input file is MS
985 #if not os.path.isdir(filename) \
986 # or not os.path.exists(filename+'/ANTENNA') \
987 # or not os.path.exists(filename+'/table.f1'):
988 if not is_ms(filename):
989 s = "File '%s' is not a Measurement set." % (filename)
990 raise IOError(s)
991 else:
992 s = "The filename should be string. "
993 raise TypeError(s)
994 # Check out put file name
995 outname=''
996 if len(outprefix) > 0: prefix=outprefix+'.'
997 else:
998 prefix=filename.rstrip('/')
999 # Now do the actual splitting.
1000 outfiles=[]
1001 tb.open(tablename=filename,nomodify=True)
1002 ant1=tb.getcol('ANTENNA1',0,-1,1)
1003 #anttab=tb.getkeyword('ANTENNA').split()[-1]
1004 anttab=tb.getkeyword('ANTENNA').lstrip('Table: ')
1005 tb.close()
1006 #tb.open(tablename=filename+'/ANTENNA',nomodify=True)
1007 tb.open(tablename=anttab,nomodify=True)
1008 nant=tb.nrows()
1009 antnames=tb.getcol('NAME',0,nant,1)
1010 tb.close()
1011 tmpname='asapmath.splitant.tmp'
1012 for antid in set(ant1):
1013 tb.open(tablename=filename,nomodify=True)
1014 tbsel=tb.query('ANTENNA1 == %s && ANTENNA2 == %s'%(antid,antid),tmpname)
1015 scan=scantable(tmpname,average=False,antenna=int(antid),getpt=getpt,freq_tolsr=freq_tolsr,)
1016 outname=prefix+antnames[antid]+'.asap'
1017 scan.save(outname,format='ASAP',overwrite=overwrite)
1018 tbsel.close()
1019 tb.close()
1020 del tbsel
1021 del scan
1022 outfiles.append(outname)
1023 os.system('rm -rf '+tmpname)
1024 del tb
1025 return outfiles
1026
1027@asaplog_post_dec
1028def _array2dOp( scan, value, mode="ADD", tsys=False, insitu=None):
1029 """
1030 This function is workaround on the basic operation of scantable
1031 with 2 dimensional float list.
1032
1033 scan: scantable operand
1034 value: float list operand
1035 mode: operation mode (ADD, SUB, MUL, DIV)
1036 tsys: if True, operate tsys as well
1037 insitu: if False, a new scantable is returned.
1038 Otherwise, the array operation is done in-sitsu.
1039 """
1040 if insitu is None: insitu = rcParams['insitu']
1041 nrow = scan.nrow()
1042 s = None
1043 from asap._asap import stmath
1044 stm = stmath()
1045 stm._setinsitu(insitu)
1046 if len( value ) == 1:
1047 s = scantable( stm._arrayop( scan, value[0], mode, tsys ) )
1048 elif len( value ) != nrow:
1049 raise ValueError( 'len(value) must be 1 or conform to scan.nrow()' )
1050 else:
1051 from asap._asap import stmath
1052 if not insitu:
1053 s = scan.copy()
1054 else:
1055 s = scan
1056 # insitu must be True as we go row by row on the same data
1057 stm._setinsitu( True )
1058 basesel = s.get_selection()
1059 # generate a new selector object based on basesel
1060 sel = selector(basesel)
1061 for irow in range( nrow ):
1062 sel.set_rows( irow )
1063 s.set_selection( sel )
1064 if len( value[irow] ) == 1:
1065 stm._unaryop( s, value[irow][0], mode, tsys )
1066 else:
1067 #stm._arrayop( s, value[irow], mode, tsys, 'channel' )
1068 stm._arrayop( s, value[irow], mode, tsys )
1069 s.set_selection(basesel)
1070 return s
Note: See TracBrowser for help on using the repository browser.