source: trunk/python/asapmath.py @ 1826

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

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