source: branches/mergetest/python/asapmath.py @ 1779

Last change on this file since 1779 was 1779, checked in by Kana Sugimoto, 14 years ago

New Development: Yes

JIRA Issue: No (test merging alma branch)

Ready for Test: Yes

Interface Changes: Yes

What Interface Changed:

Test Programs:

Put in Release Notes: No

Module(s):

Description:


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