source: trunk/python/asapmath.py@ 3003

Last change on this file since 3003 was 2952, checked in by WataruKawasaki, 10 years ago

New Development: No

JIRA Issue: Yes CAS-6598

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed: add a parameter for some functions

Test Programs: test_sdscale

Put in Release Notes:

Module(s): sd

Description: add a parameter 'skip_flaggedrow' for STMath::unaryOperate(), STMath::arrayOperate(), STMath::arrayOperateChannel(), STMath::arrayOperateRow() and their python interfaces if exist. the default value of 'skip_flaggedrow' is false so default behaviour of these functions will not change.


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