source: trunk/python/asapmath.py@ 2514

Last change on this file since 2514 was 2472, checked in by Takeshi Nakazato, 13 years ago

New Development: No

JIRA Issue: No

Ready for Test: 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...

Updated online help.


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