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

Last change on this file since 1688 was 1686, checked in by Takeshi Nakazato, 15 years ago

New Development: No

JIRA Issue: No

Ready to Release: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: List test programs

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Changed splitant() to use new scantable constructor.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 36.5 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()
[1673]874 if ( antname.find( 'APEX' ) != -1 ):
875 scal = apexcal( scantab, scannos, calmode, verify )
876 elif ( antname.find( 'ALMA' ) != -1 or antname.find( 'OSF' ) != -1 ):
[1633]877 scal = almacal( scantab, scannos, calmode, verify )
878 else:
879 scal = calps( scantab, scannos=scannos, verify=verify )
880 elif ( calmode == 'fs' or calmode == 'fsotf' ):
881 asaplog.push( 'Calibrating %s frequency-switched data.' % antname )
882 print_log()
[1673]883 if ( antname.find( 'APEX' ) != -1 ):
884 scal = apexcal( scantab, scannos, calmode, verify )
885 elif ( antname.find( 'ALMA' ) != -1 or antname.find( 'OSF' ) != -1 ):
[1633]886 scal = almacal( scantab, scannos, calmode, verify )
887 else:
888 scal = calfs( scantab, scannos=scannos, verify=verify )
889 elif ( calmode == 'otf' ):
890 asaplog.push( 'Calibrating %s On-The-Fly data.' % antname )
891 print_log()
892 scal = almacal( scantab, scannos, calmode, verify )
893 else:
894 asaplog.push( 'No calibration.' )
895 scal = scantab.copy()
896
897 return scal
898
[1673]899def apexcal( scantab, scannos=[], calmode='none', verify=False ):
[1633]900 """
[1673]901 Calibrate APEX data
[1633]902
903 Parameters:
904 scantab: scantable
905 scannos: list of scan number
906 calmode: calibration mode
[1650]907
[1633]908 verify: verify calibration
909 """
910 from asap._asap import stmath
911 stm = stmath()
912 antname = scantab.get_antennaname()
913 ssub = scantab.get_scan( scannos )
914 scal = scantable( stm.cwcal( ssub, calmode, antname ) )
915 return scal
[1650]916
[1673]917def almacal( scantab, scannos=[], calmode='none', verify=False ):
918 """
919 Calibrate ALMA data
920
921 Parameters:
922 scantab: scantable
923 scannos: list of scan number
924 calmode: calibration mode
925
926 verify: verify calibration
927 """
928 from asap._asap import stmath
929 stm = stmath()
930 ssub = scantab.get_scan( scannos )
931 scal = scantable( stm.almacal( ssub, calmode ) )
932 return scal
933
[1650]934def splitant(filename, outprefix='',overwrite=False):
935 """
936 Split Measurement set by antenna name, save data as a scantables,
937 and return a list of filename.
938 Notice this method can only be available from CASA.
939 Prameter
940 filename: the name of Measurement set to be read.
941 outprefix: the prefix of output scantable name.
942 the names of output scantable will be
943 outprefix.antenna1, outprefix.antenna2, ....
[1651]944 If not specified, outprefix = filename is assumed.
[1650]945 overwrite If the file should be overwritten if it exists.
946 The default False is to return with warning
947 without writing the output. USE WITH CARE.
948
949 """
950 # Import the table toolkit from CASA
951 try:
[1686]952 import casac
[1650]953 except ImportError:
[1686]954 if rcParams['verbose']:
955 #print "failed to load casa"
956 print_log()
957 asaplog.push("failed to load casa")
958 print_log('ERROR')
959 else: raise
960 return False
[1650]961 try:
[1686]962 tbtool = casac.homefinder.find_home_by_name('tableHome')
963 tb = tbtool.create()
964 tb2 = tbtool.create()
[1650]965 except:
[1686]966 if rcParams['verbose']:
967 #print "failed to load a table tool:\n", e
968 print_log()
969 asaplog.push("failed to load table tool")
970 print_log('ERROR')
971 else: raise
972 return False
[1650]973 # Check the input filename
974 if isinstance(filename, str):
975 import os.path
976 filename = os.path.expandvars(filename)
977 filename = os.path.expanduser(filename)
978 if not os.path.exists(filename):
979 s = "File '%s' not found." % (filename)
980 if rcParams['verbose']:
981 print_log()
982 asaplog.push(s)
983 print_log('ERROR')
984 return
985 raise IOError(s)
986 # check if input file is MS
987 if not os.path.isdir(filename) \
988 or not os.path.exists(filename+'/ANTENNA') \
989 or not os.path.exists(filename+'/table.f1'):
990 s = "File '%s' is not a Measurement set." % (filename)
991 if rcParams['verbose']:
992 print_log()
993 asaplog.push(s)
994 print_log('ERROR')
995 return
996 raise IOError(s)
997 else:
998 s = "The filename should be string. "
999 if rcParams['verbose']:
1000 print_log()
1001 asaplog.push(s)
1002 print_log('ERROR')
1003 return
1004 raise TypeError(s)
1005 # Check out put file name
1006 outname=''
1007 if len(outprefix) > 0: prefix=outprefix+'.'
[1651]1008 else:
[1686]1009 prefix=filename.rstrip('/')
[1650]1010 # Now do the actual splitting.
1011 outfiles=[]
1012 tb.open(tablename=filename+'/ANTENNA',nomodify=True)
1013 nant=tb.nrows()
1014 antnames=tb.getcol('NAME',0,nant,1)
1015 antpos=tb.getcol('POSITION',0,nant,1).transpose()
1016 tb.close()
1017 tb.open(tablename=filename,nomodify=True)
1018 ant1=tb.getcol('ANTENNA1',0,-1,1)
[1686]1019 tb.close()
[1650]1020 for antid in set(ant1):
[1686]1021 scan=scantable(filename,average=False,getpt=True,antenna=int(antid))
[1650]1022 outname=prefix+antnames[antid]+'.asap'
1023 scan.save(outname,format='ASAP',overwrite=overwrite)
[1686]1024 del scan
[1650]1025 outfiles.append(outname)
1026 del tb, tb2
1027 return outfiles
[1679]1028
1029def _array2dOp( scan, value, mode="ADD", tsys=False ):
1030 """
1031 This function is workaround on the basic operation of scantable
1032 with 2 dimensional float list.
1033
1034 scan: scantable operand
1035 value: float list operand
1036 mode: operation mode (ADD, SUB, MUL, DIV)
1037 tsys: if True, operate tsys as well
1038 """
1039 nrow = scan.nrow()
1040 s = None
[1682]1041 if len( value ) == 1:
1042 from asap._asap import stmath
1043 stm = stmath()
1044 s = scantable( stm._arrayop( scan.copy(), value[0], mode, tsys ) )
1045 del stm
1046 elif len( value ) != nrow:
1047 asaplog.push( 'len(value) must be 1 or conform to scan.nrow()' )
[1679]1048 print_log( 'ERROR' )
1049 else:
1050 from asap._asap import stmath
1051 stm = stmath()
1052 # insitu must be True
1053 stm._setinsitu( True )
1054 s = scan.copy()
1055 sel = selector()
1056 for irow in range( nrow ):
1057 sel.set_rows( irow )
1058 s.set_selection( sel )
1059 if len( value[irow] ) == 1:
1060 stm._unaryop( s, value[irow][0], mode, tsys )
1061 else:
1062 stm._arrayop( s, value[irow], mode, tsys, 'channel' )
1063 s.set_selection()
1064 sel.reset()
1065 del sel
1066 del stm
1067 return s
1068
1069
1070
Note: See TracBrowser for help on using the repository browser.