source: trunk/python/asapmath.py@ 1878

Last change on this file since 1878 was 1862, checked in by Malte Marquarding, 14 years ago

renamed print_log_dec to more explicit asaplog_post_dec

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