source: branches/alma/python/asapmath.py @ 1659

Last change on this file since 1659 was 1659, checked in by Kana Sugimoto, 15 years ago

New Development: No

JIRA Issue: Yes (CAS-1441)

Ready to Release: Yes

Interface Changes: No

What Interface Changed:

Test Programs: run (sd.)splitant with multi-point observation data.

Put in Release Notes: No

Module(s):

Description: Fixed a bug. Pointing data is now properly handled.


  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 35.4 KB
Line 
1from asap.scantable import scantable
2from asap import rcParams
3from asap import print_log
4from asap import selector
5from asap import asaplog
6from asap import asaplotgui
7
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
124def dototalpower(calon, caloff, tcalval=0.0):
125    """
126    Do calibration for CAL on,off signals.
127    Adopted from GBTIDL dototalpower
128    Parameters:
129        calon:         the 'cal on' subintegration
130        caloff:        the 'cal off' subintegration
131        tcalval:       user supplied Tcal value
132    """
133    varlist = vars()
134    from asap._asap import stmath
135    stm = stmath()
136    stm._setinsitu(False)
137    s = scantable(stm._dototalpower(calon, caloff, tcalval))
138    s._add_history("dototalpower",varlist)
139    print_log()
140    return s
141
142def dosigref(sig, ref, smooth, tsysval=0.0, tauval=0.0):
143    """
144    Calculate a quotient (sig-ref/ref * Tsys)
145    Adopted from GBTIDL dosigref
146    Parameters:
147        sig:         on source data
148        ref:         reference data
149        smooth:      width of box car smoothing for reference
150        tsysval:     user specified Tsys (scalar only)
151        tauval:      user specified Tau (required if tsysval is set)
152    """
153    varlist = vars()
154    from asap._asap import stmath
155    stm = stmath()
156    stm._setinsitu(False)
157    s = scantable(stm._dosigref(sig, ref, smooth, tsysval, tauval))
158    s._add_history("dosigref",varlist)
159    print_log()
160    return s
161
162def calps(scantab, scannos, smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
163    """
164    Calibrate GBT position switched data
165    Adopted from GBTIDL getps
166    Currently calps identify the scans as position switched data if they
167    contain '_ps' in the source name. The data must contains 'CAL' signal
168    on/off in each integration. To identify 'CAL' on state, the word, 'calon'
169    need to be present in the source name field.
170    (GBT MS data reading process to scantable automatically append these
171    id names to the source names)
172
173    Parameters:
174        scantab:       scantable
175        scannos:       list of scan numbers
176        smooth:        optional box smoothing order for the reference
177                       (default is 1 = no smoothing)
178        tsysval:       optional user specified Tsys (default is 0.0,
179                       use Tsys in the data)
180        tauval:        optional user specified Tau
181        tcalval:       optional user specified Tcal (default is 0.0,
182                       use Tcal value in the data)
183    """
184    varlist = vars()
185    # check for the appropriate data
186    s = scantab.get_scan('*_ps*')
187    if s is None:
188        msg = "The input data appear to contain no position-switch mode data."
189        if rcParams['verbose']:
190            #print msg
191            asaplog.push(msg)
192            print_log('ERROR')
193            return
194        else:
195            raise TypeError(msg)
196    ssub = s.get_scan(scannos)
197    if ssub is None:
198        msg = "No data was found with given scan numbers!"
199        if rcParams['verbose']:
200            #print msg
201            asaplog.push(msg)
202            print_log('ERROR')
203            return
204        else:
205            raise TypeError(msg)
206    ssubon = ssub.get_scan('*calon')
207    ssuboff = ssub.get_scan('*[^calon]')
208    if ssubon.nrow() != ssuboff.nrow():
209        msg = "mismatch in numbers of CAL on/off scans. Cannot calibrate. Check the scan numbers."
210        if rcParams['verbose']:
211            #print msg
212            asaplog.push(msg)
213            print_log('ERROR')
214            return
215        else:
216            raise TypeError(msg)
217    cals = dototalpower(ssubon, ssuboff, tcalval)
218    sig = cals.get_scan('*ps')
219    ref = cals.get_scan('*psr')
220    if sig.nscan() != ref.nscan():
221        msg = "mismatch in numbers of on/off scans. Cannot calibrate. Check the scan numbers."
222        if rcParams['verbose']:
223            #print msg
224            asaplog.push(msg)
225            print_log('ERROR')
226            return
227        else:
228            raise TypeError(msg)
229
230    #for user supplied Tsys
231    if tsysval>0.0:
232        if tauval<=0.0:
233            msg = "Need to supply a valid tau to use the supplied Tsys"
234            if rcParams['verbose']:
235                #print msg
236                asaplog.push(msg)
237                print_log('ERROR')
238                return
239            else:
240                raise TypeError(msg)
241        else:
242            sig.recalc_azel()
243            ref.recalc_azel()
244            #msg = "Use of user specified Tsys is not fully implemented yet."
245            #if rcParams['verbose']:
246            #    print msg
247            #    return
248            #else:
249            #    raise TypeError(msg)
250            # use get_elevation to get elevation and
251            # calculate a scaling factor using the formula
252            # -> tsys use to dosigref
253
254    #ress = dosigref(sig, ref, smooth, tsysval)
255    ress = dosigref(sig, ref, smooth, tsysval, tauval)
256    ###
257    if verify:
258        # get data
259        import numpy
260        precal={}
261        postcal=[]
262        keys=['ps','ps_calon','psr','psr_calon']
263        ifnos=list(ssub.getifnos())
264        polnos=list(ssub.getpolnos())
265        sel=selector()
266        for i in range(2):
267            ss=ssuboff.get_scan('*'+keys[2*i])
268            ll=[]
269            for j in range(len(ifnos)):
270                for k in range(len(polnos)):
271                    sel.set_ifs(ifnos[j])
272                    sel.set_polarizations(polnos[k])
273                    try:
274                        ss.set_selection(sel)
275                    except:
276                        continue
277                    ll.append(numpy.array(ss._getspectrum(0)))
278                    sel.reset()
279                    ss.set_selection()
280            precal[keys[2*i]]=ll
281            del ss
282            ss=ssubon.get_scan('*'+keys[2*i+1])
283            ll=[]
284            for j in range(len(ifnos)):
285                for k in range(len(polnos)):
286                    sel.set_ifs(ifnos[j])
287                    sel.set_polarizations(polnos[k])
288                    try:
289                        ss.set_selection(sel)
290                    except:
291                        continue
292                    ll.append(numpy.array(ss._getspectrum(0)))
293                    sel.reset()
294                    ss.set_selection()
295            precal[keys[2*i+1]]=ll
296            del ss
297        for j in range(len(ifnos)):
298            for k in range(len(polnos)):
299                sel.set_ifs(ifnos[j])
300                sel.set_polarizations(polnos[k])
301                try:
302                    ress.set_selection(sel)
303                except:
304                    continue
305                postcal.append(numpy.array(ress._getspectrum(0)))
306                sel.reset()
307                ress.set_selection()
308        del sel
309        # plot
310        print_log()
311        asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.')
312        print_log('WARN')
313        p=asaplotgui.asaplotgui()
314        #nr=min(6,len(ifnos)*len(polnos))
315        nr=len(ifnos)*len(polnos)
316        titles=[]
317        btics=[]
318        if nr<4:
319            p.set_panels(rows=nr,cols=2,nplots=2*nr,ganged=False)
320            for i in range(2*nr):
321                b=False
322                if i >= 2*nr-2:
323                    b=True
324                btics.append(b)
325        elif nr==4:
326            p.set_panels(rows=2,cols=4,nplots=8,ganged=False)
327            for i in range(2*nr):
328                b=False
329                if i >= 2*nr-4:
330                    b=True
331                btics.append(b)
332        elif nr<7:
333            p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
334            for i in range(2*nr):
335                if i >= 2*nr-4:
336                    b=True
337                btics.append(b)
338        else:
339            print_log()
340            asaplog.push('Only first 6 [if,pol] pairs are plotted.')
341            print_log('WARN')
342            nr=6
343            for i in range(2*nr):
344                b=False
345                if i >= 2*nr-4:
346                    b=True
347                btics.append(b)
348            p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
349        for i in range(nr):
350            p.subplot(2*i)
351            p.color=0
352            title='raw data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
353            titles.append(title)
354            #p.set_axes('title',title,fontsize=40)
355            ymin=1.0e100
356            ymax=-1.0e100
357            nchan=s.nchan()
358            edge=int(nchan*0.01)
359            for j in range(4):
360                spmin=min(precal[keys[j]][i][edge:nchan-edge])
361                spmax=max(precal[keys[j]][i][edge:nchan-edge])
362                ymin=min(ymin,spmin)
363                ymax=max(ymax,spmax)
364            for j in range(4):
365                if i==0:
366                    p.set_line(label=keys[j])
367                else:
368                    p.legend()
369                p.plot(precal[keys[j]][i])
370            p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
371            if not btics[2*i]:
372                p.axes.set_xticks([])
373            p.subplot(2*i+1)
374            p.color=0
375            title='cal data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
376            titles.append(title)
377            #p.set_axes('title',title)
378            p.legend()
379            ymin=postcal[i][edge:nchan-edge].min()
380            ymax=postcal[i][edge:nchan-edge].max()
381            p.plot(postcal[i])
382            p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
383            if not btics[2*i+1]:
384                p.axes.set_xticks([])
385        for i in range(2*nr):
386            p.subplot(i)
387            p.set_axes('title',titles[i],fontsize='medium')
388        x=raw_input('Accept calibration ([y]/n): ' )
389        if x.upper() == 'N':
390            p.unmap()
391            del p
392            return scabtab
393        p.unmap()
394        del p
395    ###
396    ress._add_history("calps", varlist)
397    print_log()
398    return ress
399
400def calnod(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
401    """
402    Do full (but a pair of scans at time) processing of GBT Nod data
403    calibration.
404    Adopted from  GBTIDL's getnod
405    Parameters:
406        scantab:     scantable
407        scannos:     a pair of scan numbers, or the first scan number of the pair
408        smooth:      box car smoothing order
409        tsysval:     optional user specified Tsys value
410        tauval:      optional user specified tau value (not implemented yet)
411        tcalval:     optional user specified Tcal value
412    """
413    varlist = vars()
414    from asap._asap import stmath
415    stm = stmath()
416    stm._setinsitu(False)
417
418    # check for the appropriate data
419    s = scantab.get_scan('*_nod*')
420    if s is None:
421        msg = "The input data appear to contain no Nod observing mode data."
422        if rcParams['verbose']:
423            #print msg
424            asaplog.push(msg)
425            print_log('ERROR')
426            return
427        else:
428            raise TypeError(msg)
429
430    # need check correspondance of each beam with sig-ref ...
431    # check for timestamps, scan numbers, subscan id (not available in
432    # ASAP data format...). Assume 1st scan of the pair (beam 0 - sig
433    # and beam 1 - ref...)
434    # First scan number of paired scans or list of pairs of
435    # scan numbers (has to have even number of pairs.)
436
437    #data splitting
438    scan1no = scan2no = 0
439
440    if len(scannos)==1:
441        scan1no = scannos[0]
442        scan2no = scannos[0]+1
443        pairScans = [scan1no, scan2no]
444    else:
445        #if len(scannos)>2:
446        #    msg = "calnod can only process a pair of nod scans at time."
447        #    if rcParams['verbose']:
448        #        print msg
449        #        return
450        #    else:
451        #        raise TypeError(msg)
452        #
453        #if len(scannos)==2:
454        #    scan1no = scannos[0]
455        #    scan2no = scannos[1]
456        pairScans = list(scannos)
457
458    if tsysval>0.0:
459        if tauval<=0.0:
460            msg = "Need to supply a valid tau to use the supplied Tsys"
461            if rcParams['verbose']:
462                #print msg
463                asaplog.push(msg)
464                print_log('ERROR')
465                return
466            else:
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=['nod','nod_calon']
478        ifnos=list(scantab.getifnos())
479        polnos=list(scantab.getpolnos())
480        sel=selector()
481        for i in range(2):
482            ss=scantab.get_scan('*'+keys[i])
483            ll=[]
484            ll2=[]
485            for j in range(len(ifnos)):
486                for k in range(len(polnos)):
487                    sel.set_ifs(ifnos[j])
488                    sel.set_polarizations(polnos[k])
489                    sel.set_scans(pairScans[0])
490                    try:
491                        ss.set_selection(sel)
492                    except:
493                        continue
494                    ll.append(numpy.array(ss._getspectrum(0)))
495                    sel.reset()
496                    ss.set_selection()
497                    sel.set_ifs(ifnos[j])
498                    sel.set_polarizations(polnos[k])
499                    sel.set_scans(pairScans[1])
500                    try:
501                        ss.set_selection(sel)
502                    except:
503                        ll.pop()
504                        continue
505                    ll2.append(numpy.array(ss._getspectrum(0)))
506                    sel.reset()
507                    ss.set_selection()
508            key='%s%s' %(pairScans[0],keys[i].lstrip('nod'))
509            precal[key]=ll
510            key='%s%s' %(pairScans[1],keys[i].lstrip('nod'))
511            precal[key]=ll2
512            del ss
513        keys=precal.keys()
514        for j in range(len(ifnos)):
515            for k in range(len(polnos)):
516                sel.set_ifs(ifnos[j])
517                sel.set_polarizations(polnos[k])
518                sel.set_scans(pairScans[0])
519                try:
520                    resspec.set_selection(sel)
521                except:
522                    continue
523                postcal.append(numpy.array(resspec._getspectrum(0)))
524                sel.reset()
525                resspec.set_selection()
526        del sel
527        # plot
528        print_log()
529        asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.')
530        print_log('WARN')
531        p=asaplotgui.asaplotgui()
532        #nr=min(6,len(ifnos)*len(polnos))
533        nr=len(ifnos)*len(polnos)
534        titles=[]
535        btics=[]
536        if nr<4:
537            p.set_panels(rows=nr,cols=2,nplots=2*nr,ganged=False)
538            for i in range(2*nr):
539                b=False
540                if i >= 2*nr-2:
541                    b=True
542                btics.append(b)
543        elif nr==4:
544            p.set_panels(rows=2,cols=4,nplots=8,ganged=False)
545            for i in range(2*nr):
546                b=False
547                if i >= 2*nr-4:
548                    b=True
549                btics.append(b)
550        elif nr<7:
551            p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
552            for i in range(2*nr):
553                if i >= 2*nr-4:
554                    b=True
555                btics.append(b)
556        else:
557            print_log()
558            asaplog.push('Only first 6 [if,pol] pairs are plotted.')
559            print_log('WARN')
560            nr=6
561            for i in range(2*nr):
562                b=False
563                if i >= 2*nr-4:
564                    b=True
565                btics.append(b)
566            p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False)
567        for i in range(nr):
568            p.subplot(2*i)
569            p.color=0
570            title='raw data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
571            titles.append(title)
572            #p.set_axes('title',title,fontsize=40)
573            ymin=1.0e100
574            ymax=-1.0e100
575            nchan=scantab.nchan()
576            edge=int(nchan*0.01)
577            for j in range(4):
578                spmin=min(precal[keys[j]][i][edge:nchan-edge])
579                spmax=max(precal[keys[j]][i][edge:nchan-edge])
580                ymin=min(ymin,spmin)
581                ymax=max(ymax,spmax)
582            for j in range(4):
583                if i==0:
584                    p.set_line(label=keys[j])
585                else:
586                    p.legend()
587                p.plot(precal[keys[j]][i])
588            p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
589            if not btics[2*i]:
590                p.axes.set_xticks([])
591            p.subplot(2*i+1)
592            p.color=0
593            title='cal data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)])
594            titles.append(title)
595            #p.set_axes('title',title)
596            p.legend()
597            ymin=postcal[i][edge:nchan-edge].min()
598            ymax=postcal[i][edge:nchan-edge].max()
599            p.plot(postcal[i])
600            p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
601            if not btics[2*i+1]:
602                p.axes.set_xticks([])
603        for i in range(2*nr):
604            p.subplot(i)
605            p.set_axes('title',titles[i],fontsize='medium')
606        x=raw_input('Accept calibration ([y]/n): ' )
607        if x.upper() == 'N':
608            p.unmap()
609            del p
610            return scabtab
611        p.unmap()
612        del p
613    ###
614    resspec._add_history("calnod",varlist)
615    print_log()
616    return resspec
617
618def calfs(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False):
619    """
620    Calibrate GBT frequency switched data.
621    Adopted from GBTIDL getfs.
622    Currently calfs identify the scans as frequency switched data if they
623    contain '_fs' in the source name. The data must contains 'CAL' signal
624    on/off in each integration. To identify 'CAL' on state, the word, 'calon'
625    need to be present in the source name field.
626    (GBT MS data reading via scantable automatically append these
627    id names to the source names)
628
629    Parameters:
630        scantab:       scantable
631        scannos:       list of scan numbers
632        smooth:        optional box smoothing order for the reference
633                       (default is 1 = no smoothing)
634        tsysval:       optional user specified Tsys (default is 0.0,
635                       use Tsys in the data)
636        tauval:        optional user specified Tau
637    """
638    varlist = vars()
639    from asap._asap import stmath
640    stm = stmath()
641    stm._setinsitu(False)
642
643#    check = scantab.get_scan('*_fs*')
644#    if check is None:
645#        msg = "The input data appear to contain no Nod observing mode data."
646#        if rcParams['verbose']:
647#            print msg
648#            return
649#        else:
650#            raise TypeError(msg)
651    s = scantab.get_scan(scannos)
652    del scantab
653
654    resspec = scantable(stm._dofs(s, scannos, smooth, tsysval,tauval,tcalval))
655    ###
656    if verify:
657        # get data
658        ssub = s.get_scan(scannos)
659        ssubon = ssub.get_scan('*calon')
660        ssuboff = ssub.get_scan('*[^calon]')
661        import numpy
662        precal={}
663        postcal=[]
664        keys=['fs','fs_calon','fsr','fsr_calon']
665        ifnos=list(ssub.getifnos())
666        polnos=list(ssub.getpolnos())
667        sel=selector()
668        for i in range(2):
669            ss=ssuboff.get_scan('*'+keys[2*i])
670            ll=[]
671            for j in range(len(ifnos)):
672                for k in range(len(polnos)):
673                    sel.set_ifs(ifnos[j])
674                    sel.set_polarizations(polnos[k])
675                    try:
676                        ss.set_selection(sel)
677                    except:
678                        continue
679                    ll.append(numpy.array(ss._getspectrum(0)))
680                    sel.reset()
681                    ss.set_selection()
682            precal[keys[2*i]]=ll
683            del ss
684            ss=ssubon.get_scan('*'+keys[2*i+1])
685            ll=[]
686            for j in range(len(ifnos)):
687                for k in range(len(polnos)):
688                    sel.set_ifs(ifnos[j])
689                    sel.set_polarizations(polnos[k])
690                    try:
691                        ss.set_selection(sel)
692                    except:
693                        continue
694                    ll.append(numpy.array(ss._getspectrum(0)))
695                    sel.reset()
696                    ss.set_selection()
697            precal[keys[2*i+1]]=ll
698            del ss
699        sig=resspec.get_scan('*_fs')
700        ref=resspec.get_scan('*_fsr')
701        for k in range(len(polnos)):
702            for j in range(len(ifnos)):
703                sel.set_ifs(ifnos[j])
704                sel.set_polarizations(polnos[k])
705                try:
706                    sig.set_selection(sel)
707                    postcal.append(numpy.array(sig._getspectrum(0)))
708                except:
709                    ref.set_selection(sel)
710                    postcal.append(numpy.array(ref._getspectrum(0)))
711                sel.reset()
712                resspec.set_selection()
713        del sel
714        # plot
715        print_log()
716        asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.')
717        print_log('WARN')
718        p=asaplotgui.asaplotgui()
719        #nr=min(6,len(ifnos)*len(polnos))
720        nr=len(ifnos)/2*len(polnos)
721        titles=[]
722        btics=[]
723        if nr>3:
724            print_log()
725            asaplog.push('Only first 3 [if,pol] pairs are plotted.')
726            print_log('WARN')
727            nr=3
728        p.set_panels(rows=nr,cols=3,nplots=3*nr,ganged=False)
729        for i in range(3*nr):
730            b=False
731            if i >= 3*nr-3:
732                b=True
733            btics.append(b)
734        for i in range(nr):
735            p.subplot(3*i)
736            p.color=0
737            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)])
738            titles.append(title)
739            #p.set_axes('title',title,fontsize=40)
740            ymin=1.0e100
741            ymax=-1.0e100
742            nchan=s.nchan()
743            edge=int(nchan*0.01)
744            for j in range(4):
745                spmin=min(precal[keys[j]][i][edge:nchan-edge])
746                spmax=max(precal[keys[j]][i][edge:nchan-edge])
747                ymin=min(ymin,spmin)
748                ymax=max(ymax,spmax)
749            for j in range(4):
750                if i==0:
751                    p.set_line(label=keys[j])
752                else:
753                    p.legend()
754                p.plot(precal[keys[j]][i])
755            p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
756            if not btics[3*i]:
757                p.axes.set_xticks([])
758            p.subplot(3*i+1)
759            p.color=0
760            title='sig data IF%s POL%s' % (ifnos[2*int(i/len(polnos))],polnos[i%len(polnos)])
761            titles.append(title)
762            #p.set_axes('title',title)
763            p.legend()
764            ymin=postcal[2*i][edge:nchan-edge].min()
765            ymax=postcal[2*i][edge:nchan-edge].max()
766            p.plot(postcal[2*i])
767            p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
768            if not btics[3*i+1]:
769                p.axes.set_xticks([])
770            p.subplot(3*i+2)
771            p.color=0
772            title='ref data IF%s POL%s' % (ifnos[2*int(i/len(polnos))+1],polnos[i%len(polnos)])
773            titles.append(title)
774            #p.set_axes('title',title)
775            p.legend()
776            ymin=postcal[2*i+1][edge:nchan-edge].min()
777            ymax=postcal[2*i+1][edge:nchan-edge].max()
778            p.plot(postcal[2*i+1])
779            p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax))
780            if not btics[3*i+2]:
781                p.axes.set_xticks([])
782        for i in range(3*nr):
783            p.subplot(i)
784            p.set_axes('title',titles[i],fontsize='medium')
785        x=raw_input('Accept calibration ([y]/n): ' )
786        if x.upper() == 'N':
787            p.unmap()
788            del p
789            return scabtab
790        p.unmap()
791        del p
792    ###
793    resspec._add_history("calfs",varlist)
794    print_log()
795    return resspec
796
797def simple_math(left, right, op='add', tsys=True):
798    """
799    Apply simple mathematical binary operations to two
800    scan tables,  returning the result in a new scan table.
801    The operation is applied to both the correlations and the TSys data
802    The cursor of the output scan is set to 0
803    Parameters:
804        left:          the 'left' scan
805        right:         the 'right' scan
806        op:            the operation: 'add' (default), 'sub', 'mul', 'div'
807        tsys:          if True (default) then apply the operation to Tsys
808                       as well as the data
809    """
810    #print "simple_math is deprecated use +=/* instead."
811    asaplog.push( "simple_math is deprecated use +=/* instead." )
812    print_log('WARN')
813
814def merge(*args):
815    """
816    Merge a list of scanatables, or comma-sperated scantables into one
817    scnatble.
818    Parameters:
819        A list [scan1, scan2] or scan1, scan2.
820    Example:
821        myscans = [scan1, scan2]
822        allscans = merge(myscans)
823        # or equivalent
824        sameallscans = merge(scan1, scan2)
825    """
826    varlist = vars()
827    if isinstance(args[0],list):
828        lst = tuple(args[0])
829    elif isinstance(args[0],tuple):
830        lst = args[0]
831    else:
832        lst = tuple(args)
833    varlist["args"] = "%d scantables" % len(lst)
834    # need special formatting her for history...
835    from asap._asap import stmath
836    stm = stmath()
837    for s in lst:
838        if not isinstance(s,scantable):
839            msg = "Please give a list of scantables"
840            if rcParams['verbose']:
841                #print msg
842                asaplog.push(msg)
843                print_log('ERROR')
844                return
845            else:
846                raise TypeError(msg)
847    s = scantable(stm._merge(lst))
848    s._add_history("merge", varlist)
849    print_log()
850    return s
851
852def calibrate( scantab, scannos=[], calmode='none', verify=None ):
853    """
854    Calibrate data.
855   
856    Parameters:
857        scantab:       scantable
858        scannos:       list of scan number
859        calmode:       calibration mode
860        verify:        verify calibration     
861    """
862    antname = scantab.get_antennaname()
863    if ( calmode == 'nod' ):
864        asaplog.push( 'Calibrating nod data.' )
865        print_log()
866        scal = calnod( scantab, scannos=scannos, verify=verify )
867    elif ( calmode == 'quotient' ):
868        asaplog.push( 'Calibrating using quotient.' )
869        print_log()
870        scal = scantab.auto_quotient( verify=verify )
871    elif ( calmode == 'ps' ):
872        asaplog.push( 'Calibrating %s position-switched data.' % antname )
873        print_log()
874        if ( antname.find( 'APEX' ) != -1 or antname.find( 'ALMA' ) != -1 ):
875            scal = almacal( scantab, scannos, calmode, verify )
876        else:
877            scal = calps( scantab, scannos=scannos, verify=verify )
878    elif ( calmode == 'fs' or calmode == 'fsotf' ):
879        asaplog.push( 'Calibrating %s frequency-switched data.' % antname )
880        print_log()
881        if ( antname.find( 'APEX' ) != -1 or antname.find( 'ALMA' ) != -1 ):
882            scal = almacal( scantab, scannos, calmode, verify )
883        else:
884            scal = calfs( scantab, scannos=scannos, verify=verify )
885    elif ( calmode == 'otf' ):
886        asaplog.push( 'Calibrating %s On-The-Fly data.' % antname )
887        print_log()
888        scal = almacal( scantab, scannos, calmode, verify )
889    else:
890        asaplog.push( 'No calibration.' )
891        scal = scantab.copy()
892
893    return scal
894
895def almacal( scantab, scannos=[], calmode='none', verify=False ):
896    """
897    Calibrate APEX and ALMA data
898
899    Parameters:
900        scantab:       scantable
901        scannos:       list of scan number
902        calmode:       calibration mode
903
904        verify:        verify calibration     
905    """
906    from asap._asap import stmath
907    stm = stmath()
908    antname = scantab.get_antennaname()
909    ssub = scantab.get_scan( scannos )
910    scal = scantable( stm.cwcal( ssub, calmode, antname ) )
911    return scal
912
913def splitant(filename, outprefix='',overwrite=False):
914    """
915    Split Measurement set by antenna name, save data as a scantables,
916    and return a list of filename.
917    Notice this method can only be available from CASA.
918    Prameter
919       filename:    the name of Measurement set to be read.
920       outprefix:   the prefix of output scantable name.
921                    the names of output scantable will be
922                    outprefix.antenna1, outprefix.antenna2, ....
923                    If not specified, outprefix = filename is assumed.
924       overwrite    If the file should be overwritten if it exists.
925                    The default False is to return with warning
926                    without writing the output. USE WITH CARE.
927                   
928    """
929    # Import the table toolkit from CASA
930    try:
931        import casac
932    except ImportError:
933        if rcParams['verbose']:
934            #print "failed to load casa"
935            print_log()
936            asaplog.push("failed to load casa")
937            print_log('ERROR')
938        else: raise
939        return False
940    try:
941        tbtool = casac.homefinder.find_home_by_name('tableHome')
942        tb = tbtool.create()
943        tb2 = tbtool.create()
944    except:
945        if rcParams['verbose']:
946            #print "failed to load a table tool:\n", e
947            print_log()
948            asaplog.push("failed to load table tool")
949            print_log('ERROR')
950        else: raise
951        return False
952    # Check the input filename
953    if isinstance(filename, str):
954        import os.path
955        filename = os.path.expandvars(filename)
956        filename = os.path.expanduser(filename)
957        if not os.path.exists(filename):
958            s = "File '%s' not found." % (filename)
959            if rcParams['verbose']:
960                print_log()
961                asaplog.push(s)
962                print_log('ERROR')
963                return
964            raise IOError(s)
965        # check if input file is MS
966        if not os.path.isdir(filename) \
967               or not os.path.exists(filename+'/ANTENNA') \
968               or not os.path.exists(filename+'/table.f1'):
969            s = "File '%s' is not a Measurement set." % (filename)
970            if rcParams['verbose']:
971                print_log()
972                asaplog.push(s)
973                print_log('ERROR')
974                return
975            raise IOError(s)
976    else:
977        s = "The filename should be string. "
978        if rcParams['verbose']:
979            print_log()
980            asaplog.push(s)
981            print_log('ERROR')
982            return
983        raise TypeError(s)
984    # Check out put file name
985    outname=''
986    if len(outprefix) > 0: prefix=outprefix+'.'
987    else:
988        prefix=filename
989    # Now do the actual splitting.
990    outfiles=[]
991    tmpms="temp_antsplit.ms"
992    if os.path.exists(tmpms):
993        ans=raw_input('Temporal file '+tmpms+' exists. Delete it and continue? [y/N]: ')
994        if ans.upper() == 'Y':
995            os.system('rm -rf '+tmpms)
996            asaplog.push('The file '+tmpms+' deleted.')
997        else:
998            asaplog.push('Exit without splitting.')
999            return
1000    tb.open(tablename=filename+'/ANTENNA',nomodify=True)
1001    nant=tb.nrows()
1002    antnames=tb.getcol('NAME',0,nant,1)
1003    antpos=tb.getcol('POSITION',0,nant,1).transpose()
1004    tb.close()
1005    tb.open(tablename=filename,nomodify=True)
1006    ant1=tb.getcol('ANTENNA1',0,-1,1)
1007    for antid in set(ant1):
1008        qstr="ANTENNA1 == "+str(antid)
1009        stab = tb.queryC(qstr)
1010        ctab = stab.copy(tmpms,deep=True)
1011        stab.close()
1012        ctab.close()
1013        scan=scantable(tmpms,average=False,getpt=True)
1014        outname=prefix+antnames[antid]+'.asap'
1015        scan.save(outname,format='ASAP',overwrite=overwrite)
1016        # Modify scantable header
1017        tb2.open(tablename=outname,nomodify=False)
1018        tb2.putkeyword(keyword='AntennaName',value=antnames[antid])
1019        tb2.putkeyword(keyword='AntennaPosition',value=antpos[antid])
1020        tb2.flush()
1021        tb2.close()
1022        del scan, ctab, stab
1023        outfiles.append(outname)
1024    tb.close()
1025    del tb, tb2
1026    os.system('rm -rf '+tmpms)
1027    return outfiles
Note: See TracBrowser for help on using the repository browser.