source: trunk/python/asapmath.py@ 2067

Last change on this file since 2067 was 2034, checked in by Takeshi Nakazato, 14 years ago

New Development: No

JIRA Issue: No

Ready for Test: Yes

Interface Changes: No

What Interface Changed: Please list interface changes

Test Programs: execute sd.splitant('msname')

Put in Release Notes: Yes/No

Module(s): Module Names change impacts.

Description: Describe your changes here...

Support for reference MS input that doesn't have own subtables.


  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 36.4 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
[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
[1920]161 Currently calps identify the scans as position switched data if source
162 type enum is pson or psoff. The data must contains 'CAL' signal
163 on/off in each integration. To identify 'CAL' on state, the source type
164 enum of poncal and poffcal need to be present in the source name field.
[1391]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)
[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')
[1819]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:
[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
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 ###
[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')
[1819]532 p=asaplotgui.asaplotgui()
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
576 nchan=scantab.nchan()
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':
609 p.unmap()
610 del p
611 return scabtab
612 p.unmap()
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
626 enum of foncal and foffcal need to be present in the source name field.
[1391]627 (GBT MS data reading via scantable automatically append these
628 id names to the source names)
629
630 Parameters:
631 scantab: scantable
632 scannos: list of scan numbers
633 smooth: optional box smoothing order for the reference
634 (default is 1 = no smoothing)
635 tsysval: optional user specified Tsys (default is 0.0,
636 use Tsys in the data)
637 tauval: optional user specified Tau
[1920]638 verify: Verify calibration if true
[1391]639 """
640 varlist = vars()
641 from asap._asap import stmath
[1819]642 from asap._asap import srctype
[1391]643 stm = stmath()
644 stm._setinsitu(False)
645
646# check = scantab.get_scan('*_fs*')
647# if check is None:
648# msg = "The input data appear to contain no Nod observing mode data."
[1859]649# raise TypeError(msg)
[1391]650 s = scantab.get_scan(scannos)
651 del scantab
652
653 resspec = scantable(stm._dofs(s, scannos, smooth, tsysval,tauval,tcalval))
[1819]654 ###
655 if verify:
656 # get data
657 ssub = s.get_scan(scannos)
658 #ssubon = ssub.get_scan('*calon')
659 #ssuboff = ssub.get_scan('*[^calon]')
660 sel = selector()
661 sel.set_types( [srctype.foncal,srctype.foffcal] )
662 ssub.set_selection( sel )
663 ssubon = ssub.copy()
664 ssub.set_selection()
665 sel.reset()
666 sel.set_types( [srctype.fson,srctype.fsoff] )
667 ssub.set_selection( sel )
668 ssuboff = ssub.copy()
669 ssub.set_selection()
670 sel.reset()
671 import numpy
672 precal={}
673 postcal=[]
674 keys=['fs','fs_calon','fsr','fsr_calon']
675 types=[srctype.fson,srctype.foncal,srctype.fsoff,srctype.foffcal]
676 ifnos=list(ssub.getifnos())
677 polnos=list(ssub.getpolnos())
678 for i in range(2):
679 #ss=ssuboff.get_scan('*'+keys[2*i])
680 ll=[]
681 for j in range(len(ifnos)):
682 for k in range(len(polnos)):
683 sel.set_ifs(ifnos[j])
684 sel.set_polarizations(polnos[k])
685 sel.set_types(types[2*i])
686 try:
687 #ss.set_selection(sel)
688 ssuboff.set_selection(sel)
689 except:
690 continue
691 ll.append(numpy.array(ss._getspectrum(0)))
692 sel.reset()
693 #ss.set_selection()
694 ssuboff.set_selection()
695 precal[keys[2*i]]=ll
696 #del ss
697 #ss=ssubon.get_scan('*'+keys[2*i+1])
698 ll=[]
699 for j in range(len(ifnos)):
700 for k in range(len(polnos)):
701 sel.set_ifs(ifnos[j])
702 sel.set_polarizations(polnos[k])
703 sel.set_types(types[2*i+1])
704 try:
705 #ss.set_selection(sel)
706 ssubon.set_selection(sel)
707 except:
708 continue
709 ll.append(numpy.array(ss._getspectrum(0)))
710 sel.reset()
711 #ss.set_selection()
712 ssubon.set_selection()
713 precal[keys[2*i+1]]=ll
714 #del ss
715 #sig=resspec.get_scan('*_fs')
716 #ref=resspec.get_scan('*_fsr')
717 sel.set_types( srctype.fson )
718 resspec.set_selection( sel )
719 sig=resspec.copy()
720 resspec.set_selection()
721 sel.reset()
722 sel.set_type( srctype.fsoff )
723 resspec.set_selection( sel )
724 ref=resspec.copy()
725 resspec.set_selection()
726 sel.reset()
727 for k in range(len(polnos)):
728 for j in range(len(ifnos)):
729 sel.set_ifs(ifnos[j])
730 sel.set_polarizations(polnos[k])
731 try:
732 sig.set_selection(sel)
733 postcal.append(numpy.array(sig._getspectrum(0)))
734 except:
735 ref.set_selection(sel)
736 postcal.append(numpy.array(ref._getspectrum(0)))
737 sel.reset()
738 resspec.set_selection()
739 del sel
740 # plot
[1861]741 asaplog.post()
[1819]742 asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.')
[1861]743 asaplog.post('WARN')
[1819]744 p=asaplotgui.asaplotgui()
745 #nr=min(6,len(ifnos)*len(polnos))
746 nr=len(ifnos)/2*len(polnos)
747 titles=[]
748 btics=[]
749 if nr>3:
[1861]750 asaplog.post()
[1819]751 asaplog.push('Only first 3 [if,pol] pairs are plotted.')
[1861]752 asaplog.post('WARN')
[1819]753 nr=3
754 p.set_panels(rows=nr,cols=3,nplots=3*nr,ganged=False)
755 for i in range(3*nr):
756 b=False
757 if i >= 3*nr-3:
758 b=True
759 btics.append(b)
760 for i in range(nr):
761 p.subplot(3*i)
762 p.color=0
763 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)])
764 titles.append(title)
765 #p.set_axes('title',title,fontsize=40)
766 ymin=1.0e100
767 ymax=-1.0e100
768 nchan=s.nchan()
769 edge=int(nchan*0.01)
770 for j in range(4):
771 spmin=min(precal[keys[j]][i][edge:nchan-edge])
772 spmax=max(precal[keys[j]][i][edge:nchan-edge])
773 ymin=min(ymin,spmin)
774 ymax=max(ymax,spmax)
775 for j in range(4):
776 if i==0:
777 p.set_line(label=keys[j])
778 else:
779 p.legend()
780 p.plot(precal[keys[j]][i])
781 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
782 if not btics[3*i]:
783 p.axes.set_xticks([])
784 p.subplot(3*i+1)
785 p.color=0
786 title='sig data IF%s POL%s' % (ifnos[2*int(i/len(polnos))],polnos[i%len(polnos)])
787 titles.append(title)
788 #p.set_axes('title',title)
789 p.legend()
790 ymin=postcal[2*i][edge:nchan-edge].min()
791 ymax=postcal[2*i][edge:nchan-edge].max()
792 p.plot(postcal[2*i])
793 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
794 if not btics[3*i+1]:
795 p.axes.set_xticks([])
796 p.subplot(3*i+2)
797 p.color=0
798 title='ref data IF%s POL%s' % (ifnos[2*int(i/len(polnos))+1],polnos[i%len(polnos)])
799 titles.append(title)
800 #p.set_axes('title',title)
801 p.legend()
802 ymin=postcal[2*i+1][edge:nchan-edge].min()
803 ymax=postcal[2*i+1][edge:nchan-edge].max()
804 p.plot(postcal[2*i+1])
805 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
806 if not btics[3*i+2]:
807 p.axes.set_xticks([])
808 for i in range(3*nr):
809 p.subplot(i)
810 p.set_axes('title',titles[i],fontsize='medium')
811 x=raw_input('Accept calibration ([y]/n): ' )
812 if x.upper() == 'N':
813 p.unmap()
814 del p
815 return scabtab
816 p.unmap()
817 del p
818 ###
[1391]819 resspec._add_history("calfs",varlist)
820 return resspec
821
[1862]822@asaplog_post_dec
[918]823def merge(*args):
[945]824 """
[1362]825 Merge a list of scanatables, or comma-sperated scantables into one
826 scnatble.
827 Parameters:
828 A list [scan1, scan2] or scan1, scan2.
829 Example:
830 myscans = [scan1, scan2]
[1589]831 allscans = merge(myscans)
832 # or equivalent
833 sameallscans = merge(scan1, scan2)
[945]834 """
[918]835 varlist = vars()
836 if isinstance(args[0],list):
837 lst = tuple(args[0])
838 elif isinstance(args[0],tuple):
839 lst = args[0]
840 else:
841 lst = tuple(args)
842 varlist["args"] = "%d scantables" % len(lst)
843 # need special formatting her for history...
844 from asap._asap import stmath
845 stm = stmath()
846 for s in lst:
847 if not isinstance(s,scantable):
848 msg = "Please give a list of scantables"
[1859]849 raise TypeError(msg)
[918]850 s = scantable(stm._merge(lst))
851 s._add_history("merge", varlist)
852 return s
[1819]853
[1862]854@asaplog_post_dec
[1819]855def calibrate( scantab, scannos=[], calmode='none', verify=None ):
856 """
857 Calibrate data.
[1826]858
[1819]859 Parameters:
860 scantab: scantable
861 scannos: list of scan number
862 calmode: calibration mode
[1826]863 verify: verify calibration
[1819]864 """
865 antname = scantab.get_antennaname()
866 if ( calmode == 'nod' ):
867 asaplog.push( 'Calibrating nod data.' )
868 scal = calnod( scantab, scannos=scannos, verify=verify )
869 elif ( calmode == 'quotient' ):
870 asaplog.push( 'Calibrating using quotient.' )
871 scal = scantab.auto_quotient( verify=verify )
872 elif ( calmode == 'ps' ):
873 asaplog.push( 'Calibrating %s position-switched data.' % antname )
874 if ( antname.find( 'APEX' ) != -1 ):
875 scal = apexcal( scantab, scannos, calmode, verify )
876 elif ( antname.find( 'ALMA' ) != -1 or antname.find( 'OSF' ) != -1 ):
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 if ( antname.find( 'APEX' ) != -1 ):
883 scal = apexcal( scantab, scannos, calmode, verify )
884 elif ( antname.find( 'ALMA' ) != -1 or antname.find( 'OSF' ) != -1 ):
885 scal = almacal( scantab, scannos, calmode, verify )
886 else:
887 scal = calfs( scantab, scannos=scannos, verify=verify )
888 elif ( calmode == 'otf' ):
889 asaplog.push( 'Calibrating %s On-The-Fly data.' % antname )
890 scal = almacal( scantab, scannos, calmode, verify )
891 else:
892 asaplog.push( 'No calibration.' )
893 scal = scantab.copy()
894
[1826]895 return scal
[1819]896
897def apexcal( scantab, scannos=[], calmode='none', verify=False ):
898 """
899 Calibrate APEX data
900
901 Parameters:
902 scantab: scantable
903 scannos: list of scan number
904 calmode: calibration mode
905
[1826]906 verify: verify calibration
[1819]907 """
908 from asap._asap import stmath
909 stm = stmath()
910 antname = scantab.get_antennaname()
911 ssub = scantab.get_scan( scannos )
912 scal = scantable( stm.cwcal( ssub, calmode, antname ) )
913 return scal
914
915def almacal( scantab, scannos=[], calmode='none', verify=False ):
916 """
917 Calibrate ALMA data
918
919 Parameters:
920 scantab: scantable
921 scannos: list of scan number
922 calmode: calibration mode
923
[1826]924 verify: verify calibration
[1819]925 """
926 from asap._asap import stmath
927 stm = stmath()
928 ssub = scantab.get_scan( scannos )
929 scal = scantable( stm.almacal( ssub, calmode ) )
930 return scal
931
[1862]932@asaplog_post_dec
[1819]933def splitant(filename, outprefix='',overwrite=False):
934 """
935 Split Measurement set by antenna name, save data as a scantables,
936 and return a list of filename.
[1826]937 Notice this method can only be available from CASA.
[1819]938 Prameter
[1826]939 filename: the name of Measurement set to be read.
[1819]940 outprefix: the prefix of output scantable name.
941 the names of output scantable will be
942 outprefix.antenna1, outprefix.antenna2, ....
943 If not specified, outprefix = filename is assumed.
944 overwrite If the file should be overwritten if it exists.
945 The default False is to return with warning
946 without writing the output. USE WITH CARE.
[1826]947
[1819]948 """
949 # Import the table toolkit from CASA
[1859]950 import casac
[1918]951 from asap.scantable import is_ms
[1859]952 tbtool = casac.homefinder.find_home_by_name('tableHome')
953 tb = tbtool.create()
[1819]954 # Check the input filename
955 if isinstance(filename, str):
956 import os.path
957 filename = os.path.expandvars(filename)
958 filename = os.path.expanduser(filename)
959 if not os.path.exists(filename):
960 s = "File '%s' not found." % (filename)
961 raise IOError(s)
962 # check if input file is MS
[1883]963 #if not os.path.isdir(filename) \
964 # or not os.path.exists(filename+'/ANTENNA') \
965 # or not os.path.exists(filename+'/table.f1'):
966 if not is_ms(filename):
[1819]967 s = "File '%s' is not a Measurement set." % (filename)
968 raise IOError(s)
969 else:
970 s = "The filename should be string. "
971 raise TypeError(s)
972 # Check out put file name
973 outname=''
974 if len(outprefix) > 0: prefix=outprefix+'.'
975 else:
976 prefix=filename.rstrip('/')
977 # Now do the actual splitting.
978 outfiles=[]
[2034]979 tb.open(tablename=filename,nomodify=True)
980 ant1=tb.getcol('ANTENNA1',0,-1,1)
981 anttab=tb.getkeyword('ANTENNA').split()[-1]
982 tb.close()
983 #tb.open(tablename=filename+'/ANTENNA',nomodify=True)
984 tb.open(tablename=anttab,nomodify=True)
[1819]985 nant=tb.nrows()
986 antnames=tb.getcol('NAME',0,nant,1)
987 tb.close()
[1880]988 tmpname='asapmath.splitant.tmp'
[1819]989 for antid in set(ant1):
[1883]990 tb.open(tablename=filename,nomodify=True)
991 tbsel=tb.query('ANTENNA1 == %s && ANTENNA2 == %s'%(antid,antid),tmpname)
[1918]992 scan=scantable(tmpname,average=False,getpt=True,antenna=int(antid))
993 outname=prefix+antnames[antid]+'.asap'
994 scan.save(outname,format='ASAP',overwrite=overwrite)
[1880]995 tbsel.close()
[1883]996 tb.close()
[1880]997 del tbsel
[1819]998 del scan
999 outfiles.append(outname)
[1880]1000 os.system('rm -rf '+tmpname)
1001 del tb
[1819]1002 return outfiles
1003
[1862]1004@asaplog_post_dec
[1819]1005def _array2dOp( scan, value, mode="ADD", tsys=False ):
1006 """
1007 This function is workaround on the basic operation of scantable
1008 with 2 dimensional float list.
1009
1010 scan: scantable operand
1011 value: float list operand
1012 mode: operation mode (ADD, SUB, MUL, DIV)
[1826]1013 tsys: if True, operate tsys as well
[1819]1014 """
1015 nrow = scan.nrow()
1016 s = None
1017 if len( value ) == 1:
1018 from asap._asap import stmath
1019 stm = stmath()
1020 s = scantable( stm._arrayop( scan.copy(), value[0], mode, tsys ) )
1021 del stm
1022 elif len( value ) != nrow:
[1859]1023 raise ValueError( 'len(value) must be 1 or conform to scan.nrow()' )
[1819]1024 else:
1025 from asap._asap import stmath
1026 stm = stmath()
[1826]1027 # insitu must be True
[1819]1028 stm._setinsitu( True )
1029 s = scan.copy()
1030 sel = selector()
1031 for irow in range( nrow ):
1032 sel.set_rows( irow )
1033 s.set_selection( sel )
1034 if len( value[irow] ) == 1:
1035 stm._unaryop( s, value[irow][0], mode, tsys )
1036 else:
1037 stm._arrayop( s, value[irow], mode, tsys, 'channel' )
1038 s.set_selection()
1039 sel.reset()
1040 del sel
1041 del stm
1042 return s
Note: See TracBrowser for help on using the repository browser.