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

Last change on this file since 1666 was 1659, checked in by Kana Sugimoto, 15 years ago

New Development: No

JIRA Issue: Yes (CAS-1441)

Ready to Release: Yes

Interface Changes: No

What Interface Changed:

Test Programs: run (sd.)splitant with multi-point observation data.

Put in Release Notes: No

Module(s):

Description: Fixed a bug. Pointing data is now properly handled.


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