|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 from asap.scantable import scantable 00002 from asap.parameters import rcParams 00003 from asap.logging import asaplog, asaplog_post_dec 00004 from asap.selector import selector 00005 from asap.asapplotter import new_asaplot 00006 from matplotlib import rc as rcp 00007 00008 @asaplog_post_dec 00009 def average_time(*args, **kwargs): 00010 """ 00011 Return the (time) average of a scan or list of scans. [in channels only] 00012 The cursor of the output scan is set to 0 00013 Parameters: 00014 one scan or comma separated scans or a list of scans 00015 mask: an optional mask (only used for 'var' and 'tsys' weighting) 00016 scanav: True averages each scan separately. 00017 False (default) averages all scans together, 00018 weight: Weighting scheme. 00019 'none' (mean no weight) 00020 'var' (1/var(spec) weighted) 00021 'tsys' (1/Tsys**2 weighted) 00022 'tint' (integration time weighted) 00023 'tintsys' (Tint/Tsys**2) 00024 'median' ( median averaging) 00025 align: align the spectra in velocity before averaging. It takes 00026 the time of the first spectrum in the first scantable 00027 as reference time. 00028 compel: True forces to average overwrapped IFs. 00029 Example: 00030 # return a time averaged scan from scana and scanb 00031 # without using a mask 00032 scanav = average_time(scana,scanb) 00033 # or equivalent 00034 # scanav = average_time([scana, scanb]) 00035 # return the (time) averaged scan, i.e. the average of 00036 # all correlator cycles 00037 scanav = average_time(scan, scanav=True) 00038 """ 00039 scanav = False 00040 if kwargs.has_key('scanav'): 00041 scanav = kwargs.get('scanav') 00042 weight = 'tint' 00043 if kwargs.has_key('weight'): 00044 weight = kwargs.get('weight') 00045 mask = () 00046 if kwargs.has_key('mask'): 00047 mask = kwargs.get('mask') 00048 align = False 00049 if kwargs.has_key('align'): 00050 align = kwargs.get('align') 00051 compel = False 00052 if kwargs.has_key('compel'): 00053 compel = kwargs.get('compel') 00054 varlist = vars() 00055 if isinstance(args[0],list): 00056 lst = args[0] 00057 elif isinstance(args[0],tuple): 00058 lst = list(args[0]) 00059 else: 00060 lst = list(args) 00061 00062 del varlist["kwargs"] 00063 varlist["args"] = "%d scantables" % len(lst) 00064 # need special formatting here for history... 00065 00066 from asap._asap import stmath 00067 stm = stmath() 00068 for s in lst: 00069 if not isinstance(s,scantable): 00070 msg = "Please give a list of scantables" 00071 raise TypeError(msg) 00072 if scanav: scanav = "SCAN" 00073 else: scanav = "NONE" 00074 alignedlst = [] 00075 if align: 00076 refepoch = lst[0].get_time(0) 00077 for scan in lst: 00078 alignedlst.append(scan.freq_align(refepoch,insitu=False)) 00079 else: 00080 alignedlst = lst 00081 if weight.upper() == 'MEDIAN': 00082 # median doesn't support list of scantables - merge first 00083 merged = None 00084 if len(alignedlst) > 1: 00085 merged = merge(alignedlst) 00086 else: 00087 merged = alignedlst[0] 00088 s = scantable(stm._averagechannel(merged, 'MEDIAN', scanav)) 00089 del merged 00090 else: 00091 #s = scantable(stm._average(alignedlst, mask, weight.upper(), scanav)) 00092 s = scantable(stm._new_average(alignedlst, compel, mask, weight.upper(), scanav)) 00093 s._add_history("average_time",varlist) 00094 00095 return s 00096 00097 @asaplog_post_dec 00098 def quotient(source, reference, preserve=True): 00099 """ 00100 Return the quotient of a 'source' (signal) scan and a 'reference' scan. 00101 The reference can have just one scan, even if the signal has many. Otherwise 00102 they must have the same number of scans. 00103 The cursor of the output scan is set to 0 00104 Parameters: 00105 source: the 'on' scan 00106 reference: the 'off' scan 00107 preserve: you can preserve (default) the continuum or 00108 remove it. The equations used are 00109 preserve: Output = Toff * (on/off) - Toff 00110 remove: Output = Toff * (on/off) - Ton 00111 """ 00112 varlist = vars() 00113 from asap._asap import stmath 00114 stm = stmath() 00115 stm._setinsitu(False) 00116 s = scantable(stm._quotient(source, reference, preserve)) 00117 s._add_history("quotient",varlist) 00118 return s 00119 00120 @asaplog_post_dec 00121 def dototalpower(calon, caloff, tcalval=0.0): 00122 """ 00123 Do calibration for CAL on,off signals. 00124 Adopted from GBTIDL dototalpower 00125 Parameters: 00126 calon: the 'cal on' subintegration 00127 caloff: the 'cal off' subintegration 00128 tcalval: user supplied Tcal value 00129 """ 00130 varlist = vars() 00131 from asap._asap import stmath 00132 stm = stmath() 00133 stm._setinsitu(False) 00134 s = scantable(stm._dototalpower(calon, caloff, tcalval)) 00135 s._add_history("dototalpower",varlist) 00136 return s 00137 00138 @asaplog_post_dec 00139 def dosigref(sig, ref, smooth, tsysval=0.0, tauval=0.0): 00140 """ 00141 Calculate a quotient (sig-ref/ref * Tsys) 00142 Adopted from GBTIDL dosigref 00143 Parameters: 00144 sig: on source data 00145 ref: reference data 00146 smooth: width of box car smoothing for reference 00147 tsysval: user specified Tsys (scalar only) 00148 tauval: user specified Tau (required if tsysval is set) 00149 """ 00150 varlist = vars() 00151 from asap._asap import stmath 00152 stm = stmath() 00153 stm._setinsitu(False) 00154 s = scantable(stm._dosigref(sig, ref, smooth, tsysval, tauval)) 00155 s._add_history("dosigref",varlist) 00156 return s 00157 00158 @asaplog_post_dec 00159 def calps(scantab, scannos, smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False): 00160 """ 00161 Calibrate GBT position switched data 00162 Adopted from GBTIDL getps 00163 Currently calps identify the scans as position switched data if source 00164 type enum is pson or psoff. The data must contains 'CAL' signal 00165 on/off in each integration. To identify 'CAL' on state, the source type 00166 enum of poncal and poffcal need to be present. 00167 00168 Parameters: 00169 scantab: scantable 00170 scannos: list of scan numbers 00171 smooth: optional box smoothing order for the reference 00172 (default is 1 = no smoothing) 00173 tsysval: optional user specified Tsys (default is 0.0, 00174 use Tsys in the data) 00175 tauval: optional user specified Tau 00176 tcalval: optional user specified Tcal (default is 0.0, 00177 use Tcal value in the data) 00178 verify: Verify calibration if true 00179 """ 00180 varlist = vars() 00181 # check for the appropriate data 00182 ## s = scantab.get_scan('*_ps*') 00183 ## if s is None: 00184 ## msg = "The input data appear to contain no position-switch mode data." 00185 ## raise TypeError(msg) 00186 s = scantab.copy() 00187 from asap._asap import srctype 00188 sel = selector() 00189 sel.set_types( srctype.pson ) 00190 try: 00191 scantab.set_selection( sel ) 00192 except Exception, e: 00193 msg = "The input data appear to contain no position-switch mode data." 00194 raise TypeError(msg) 00195 s.set_selection() 00196 sel.reset() 00197 ssub = s.get_scan(scannos) 00198 if ssub is None: 00199 msg = "No data was found with given scan numbers!" 00200 raise TypeError(msg) 00201 #ssubon = ssub.get_scan('*calon') 00202 #ssuboff = ssub.get_scan('*[^calon]') 00203 sel.set_types( [srctype.poncal,srctype.poffcal] ) 00204 ssub.set_selection( sel ) 00205 ssubon = ssub.copy() 00206 ssub.set_selection() 00207 sel.reset() 00208 sel.set_types( [srctype.pson,srctype.psoff] ) 00209 ssub.set_selection( sel ) 00210 ssuboff = ssub.copy() 00211 ssub.set_selection() 00212 sel.reset() 00213 if ssubon.nrow() != ssuboff.nrow(): 00214 msg = "mismatch in numbers of CAL on/off scans. Cannot calibrate. Check the scan numbers." 00215 raise TypeError(msg) 00216 cals = dototalpower(ssubon, ssuboff, tcalval) 00217 #sig = cals.get_scan('*ps') 00218 #ref = cals.get_scan('*psr') 00219 sel.set_types( srctype.pson ) 00220 cals.set_selection( sel ) 00221 sig = cals.copy() 00222 cals.set_selection() 00223 sel.reset() 00224 sel.set_types( srctype.psoff ) 00225 cals.set_selection( sel ) 00226 ref = cals.copy() 00227 cals.set_selection() 00228 sel.reset() 00229 if sig.nscan() != ref.nscan(): 00230 msg = "mismatch in numbers of on/off scans. Cannot calibrate. Check the scan numbers." 00231 raise TypeError(msg) 00232 00233 #for user supplied Tsys 00234 if tsysval>0.0: 00235 if tauval<=0.0: 00236 msg = "Need to supply a valid tau to use the supplied Tsys" 00237 raise TypeError(msg) 00238 else: 00239 sig.recalc_azel() 00240 ref.recalc_azel() 00241 #msg = "Use of user specified Tsys is not fully implemented yet." 00242 #raise TypeError(msg) 00243 # use get_elevation to get elevation and 00244 # calculate a scaling factor using the formula 00245 # -> tsys use to dosigref 00246 00247 #ress = dosigref(sig, ref, smooth, tsysval) 00248 ress = dosigref(sig, ref, smooth, tsysval, tauval) 00249 ### 00250 if verify: 00251 # get data 00252 import numpy 00253 precal={} 00254 postcal=[] 00255 keys=['ps','ps_calon','psr','psr_calon'] 00256 types=[srctype.pson,srctype.poncal,srctype.psoff,srctype.poffcal] 00257 ifnos=list(ssub.getifnos()) 00258 polnos=list(ssub.getpolnos()) 00259 sel=selector() 00260 for i in range(2): 00261 #ss=ssuboff.get_scan('*'+keys[2*i]) 00262 ll=[] 00263 for j in range(len(ifnos)): 00264 for k in range(len(polnos)): 00265 sel.set_ifs(ifnos[j]) 00266 sel.set_polarizations(polnos[k]) 00267 sel.set_types(types[2*i]) 00268 try: 00269 #ss.set_selection(sel) 00270 ssuboff.set_selection(sel) 00271 except: 00272 continue 00273 #ll.append(numpy.array(ss._getspectrum(0))) 00274 ll.append(numpy.array(ssuboff._getspectrum(0))) 00275 sel.reset() 00276 ssuboff.set_selection() 00277 precal[keys[2*i]]=ll 00278 #del ss 00279 #ss=ssubon.get_scan('*'+keys[2*i+1]) 00280 ll=[] 00281 for j in range(len(ifnos)): 00282 for k in range(len(polnos)): 00283 sel.set_ifs(ifnos[j]) 00284 sel.set_polarizations(polnos[k]) 00285 sel.set_types(types[2*i+1]) 00286 try: 00287 #ss.set_selection(sel) 00288 ssubon.set_selection(sel) 00289 except: 00290 continue 00291 #ll.append(numpy.array(ss._getspectrum(0))) 00292 ll.append(numpy.array(ssubon._getspectrum(0))) 00293 sel.reset() 00294 ssubon.set_selection() 00295 precal[keys[2*i+1]]=ll 00296 #del ss 00297 for j in range(len(ifnos)): 00298 for k in range(len(polnos)): 00299 sel.set_ifs(ifnos[j]) 00300 sel.set_polarizations(polnos[k]) 00301 try: 00302 ress.set_selection(sel) 00303 except: 00304 continue 00305 postcal.append(numpy.array(ress._getspectrum(0))) 00306 sel.reset() 00307 ress.set_selection() 00308 del sel 00309 # plot 00310 asaplog.post() 00311 asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.') 00312 asaplog.post('WARN') 00313 p=new_asaplot() 00314 rcp('lines', linewidth=1) 00315 #nr=min(6,len(ifnos)*len(polnos)) 00316 nr=len(ifnos)*len(polnos) 00317 titles=[] 00318 btics=[] 00319 if nr<4: 00320 p.set_panels(rows=nr,cols=2,nplots=2*nr,ganged=False) 00321 for i in range(2*nr): 00322 b=False 00323 if i >= 2*nr-2: 00324 b=True 00325 btics.append(b) 00326 elif nr==4: 00327 p.set_panels(rows=2,cols=4,nplots=8,ganged=False) 00328 for i in range(2*nr): 00329 b=False 00330 if i >= 2*nr-4: 00331 b=True 00332 btics.append(b) 00333 elif nr<7: 00334 p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False) 00335 for i in range(2*nr): 00336 if i >= 2*nr-4: 00337 b=True 00338 btics.append(b) 00339 else: 00340 asaplog.post() 00341 asaplog.push('Only first 6 [if,pol] pairs are plotted.') 00342 asaplog.post('WARN') 00343 nr=6 00344 for i in range(2*nr): 00345 b=False 00346 if i >= 2*nr-4: 00347 b=True 00348 btics.append(b) 00349 p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False) 00350 for i in range(nr): 00351 p.subplot(2*i) 00352 p.color=0 00353 title='raw data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)]) 00354 titles.append(title) 00355 #p.set_axes('title',title,fontsize=40) 00356 ymin=1.0e100 00357 ymax=-1.0e100 00358 nchan=s.nchan(ifnos[int(i/len(polnos))]) 00359 edge=int(nchan*0.01) 00360 for j in range(4): 00361 spmin=min(precal[keys[j]][i][edge:nchan-edge]) 00362 spmax=max(precal[keys[j]][i][edge:nchan-edge]) 00363 ymin=min(ymin,spmin) 00364 ymax=max(ymax,spmax) 00365 for j in range(4): 00366 if i==0: 00367 p.set_line(label=keys[j]) 00368 else: 00369 p.legend() 00370 p.plot(precal[keys[j]][i]) 00371 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax)) 00372 if not btics[2*i]: 00373 p.axes.set_xticks([]) 00374 p.subplot(2*i+1) 00375 p.color=0 00376 title='cal data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)]) 00377 titles.append(title) 00378 #p.set_axes('title',title) 00379 p.legend() 00380 ymin=postcal[i][edge:nchan-edge].min() 00381 ymax=postcal[i][edge:nchan-edge].max() 00382 p.plot(postcal[i]) 00383 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax)) 00384 if not btics[2*i+1]: 00385 p.axes.set_xticks([]) 00386 for i in range(2*nr): 00387 p.subplot(i) 00388 p.set_axes('title',titles[i],fontsize='medium') 00389 x=raw_input('Accept calibration ([y]/n): ' ) 00390 if x.upper() == 'N': 00391 p.quit() 00392 del p 00393 return scabtab 00394 p.quit() 00395 del p 00396 ### 00397 ress._add_history("calps", varlist) 00398 return ress 00399 00400 @asaplog_post_dec 00401 def calnod(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False): 00402 """ 00403 Do full (but a pair of scans at time) processing of GBT Nod data 00404 calibration. 00405 Adopted from GBTIDL's getnod 00406 Parameters: 00407 scantab: scantable 00408 scannos: a pair of scan numbers, or the first scan number of the pair 00409 smooth: box car smoothing order 00410 tsysval: optional user specified Tsys value 00411 tauval: optional user specified tau value (not implemented yet) 00412 tcalval: optional user specified Tcal value 00413 verify: Verify calibration if true 00414 """ 00415 varlist = vars() 00416 from asap._asap import stmath 00417 from asap._asap import srctype 00418 stm = stmath() 00419 stm._setinsitu(False) 00420 00421 # check for the appropriate data 00422 ## s = scantab.get_scan('*_nod*') 00423 ## if s is None: 00424 ## msg = "The input data appear to contain no Nod observing mode data." 00425 ## raise TypeError(msg) 00426 s = scantab.copy() 00427 sel = selector() 00428 sel.set_types( srctype.nod ) 00429 try: 00430 s.set_selection( sel ) 00431 except Exception, e: 00432 msg = "The input data appear to contain no Nod observing mode data." 00433 raise TypeError(msg) 00434 sel.reset() 00435 del sel 00436 del s 00437 00438 # need check correspondance of each beam with sig-ref ... 00439 # check for timestamps, scan numbers, subscan id (not available in 00440 # ASAP data format...). Assume 1st scan of the pair (beam 0 - sig 00441 # and beam 1 - ref...) 00442 # First scan number of paired scans or list of pairs of 00443 # scan numbers (has to have even number of pairs.) 00444 00445 #data splitting 00446 scan1no = scan2no = 0 00447 00448 if len(scannos)==1: 00449 scan1no = scannos[0] 00450 scan2no = scannos[0]+1 00451 pairScans = [scan1no, scan2no] 00452 else: 00453 #if len(scannos)>2: 00454 # msg = "calnod can only process a pair of nod scans at time." 00455 # raise TypeError(msg) 00456 # 00457 #if len(scannos)==2: 00458 # scan1no = scannos[0] 00459 # scan2no = scannos[1] 00460 pairScans = list(scannos) 00461 00462 if tsysval>0.0: 00463 if tauval<=0.0: 00464 msg = "Need to supply a valid tau to use the supplied Tsys" 00465 raise TypeError(msg) 00466 else: 00467 scantab.recalc_azel() 00468 resspec = scantable(stm._donod(scantab, pairScans, smooth, tsysval,tauval,tcalval)) 00469 ### 00470 if verify: 00471 # get data 00472 import numpy 00473 precal={} 00474 postcal=[] 00475 keys=['','_calon'] 00476 types=[srctype.nod,srctype.nodcal] 00477 ifnos=list(scantab.getifnos()) 00478 polnos=list(scantab.getpolnos()) 00479 sel=selector() 00480 ss = scantab.copy() 00481 for i in range(2): 00482 #ss=scantab.get_scan('*'+keys[i]) 00483 ll=[] 00484 ll2=[] 00485 for j in range(len(ifnos)): 00486 for k in range(len(polnos)): 00487 sel.set_ifs(ifnos[j]) 00488 sel.set_polarizations(polnos[k]) 00489 sel.set_scans(pairScans[0]) 00490 sel.set_types(types[i]) 00491 try: 00492 ss.set_selection(sel) 00493 except: 00494 continue 00495 ll.append(numpy.array(ss._getspectrum(0))) 00496 sel.reset() 00497 ss.set_selection() 00498 sel.set_ifs(ifnos[j]) 00499 sel.set_polarizations(polnos[k]) 00500 sel.set_scans(pairScans[1]) 00501 sel.set_types(types[i]) 00502 try: 00503 ss.set_selection(sel) 00504 except: 00505 ll.pop() 00506 continue 00507 ll2.append(numpy.array(ss._getspectrum(0))) 00508 sel.reset() 00509 ss.set_selection() 00510 key='%s%s' %(pairScans[0],keys[i]) 00511 precal[key]=ll 00512 key='%s%s' %(pairScans[1],keys[i]) 00513 precal[key]=ll2 00514 #del ss 00515 keys=precal.keys() 00516 for j in range(len(ifnos)): 00517 for k in range(len(polnos)): 00518 sel.set_ifs(ifnos[j]) 00519 sel.set_polarizations(polnos[k]) 00520 sel.set_scans(pairScans[0]) 00521 try: 00522 resspec.set_selection(sel) 00523 except: 00524 continue 00525 postcal.append(numpy.array(resspec._getspectrum(0))) 00526 sel.reset() 00527 resspec.set_selection() 00528 del sel 00529 # plot 00530 asaplog.post() 00531 asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.') 00532 asaplog.post('WARN') 00533 p=new_asaplot() 00534 rcp('lines', linewidth=1) 00535 #nr=min(6,len(ifnos)*len(polnos)) 00536 nr=len(ifnos)*len(polnos) 00537 titles=[] 00538 btics=[] 00539 if nr<4: 00540 p.set_panels(rows=nr,cols=2,nplots=2*nr,ganged=False) 00541 for i in range(2*nr): 00542 b=False 00543 if i >= 2*nr-2: 00544 b=True 00545 btics.append(b) 00546 elif nr==4: 00547 p.set_panels(rows=2,cols=4,nplots=8,ganged=False) 00548 for i in range(2*nr): 00549 b=False 00550 if i >= 2*nr-4: 00551 b=True 00552 btics.append(b) 00553 elif nr<7: 00554 p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False) 00555 for i in range(2*nr): 00556 if i >= 2*nr-4: 00557 b=True 00558 btics.append(b) 00559 else: 00560 asaplog.post() 00561 asaplog.push('Only first 6 [if,pol] pairs are plotted.') 00562 asaplog.post('WARN') 00563 nr=6 00564 for i in range(2*nr): 00565 b=False 00566 if i >= 2*nr-4: 00567 b=True 00568 btics.append(b) 00569 p.set_panels(rows=3,cols=4,nplots=2*nr,ganged=False) 00570 for i in range(nr): 00571 p.subplot(2*i) 00572 p.color=0 00573 title='raw data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)]) 00574 titles.append(title) 00575 #p.set_axes('title',title,fontsize=40) 00576 ymin=1.0e100 00577 ymax=-1.0e100 00578 nchan=scantab.nchan(ifnos[int(i/len(polnos))]) 00579 edge=int(nchan*0.01) 00580 for j in range(4): 00581 spmin=min(precal[keys[j]][i][edge:nchan-edge]) 00582 spmax=max(precal[keys[j]][i][edge:nchan-edge]) 00583 ymin=min(ymin,spmin) 00584 ymax=max(ymax,spmax) 00585 for j in range(4): 00586 if i==0: 00587 p.set_line(label=keys[j]) 00588 else: 00589 p.legend() 00590 p.plot(precal[keys[j]][i]) 00591 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax)) 00592 if not btics[2*i]: 00593 p.axes.set_xticks([]) 00594 p.subplot(2*i+1) 00595 p.color=0 00596 title='cal data IF%s POL%s' % (ifnos[int(i/len(polnos))],polnos[i%len(polnos)]) 00597 titles.append(title) 00598 #p.set_axes('title',title) 00599 p.legend() 00600 ymin=postcal[i][edge:nchan-edge].min() 00601 ymax=postcal[i][edge:nchan-edge].max() 00602 p.plot(postcal[i]) 00603 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax)) 00604 if not btics[2*i+1]: 00605 p.axes.set_xticks([]) 00606 for i in range(2*nr): 00607 p.subplot(i) 00608 p.set_axes('title',titles[i],fontsize='medium') 00609 x=raw_input('Accept calibration ([y]/n): ' ) 00610 if x.upper() == 'N': 00611 p.quit() 00612 del p 00613 return scabtab 00614 p.quit() 00615 del p 00616 ### 00617 resspec._add_history("calnod",varlist) 00618 return resspec 00619 00620 @asaplog_post_dec 00621 def calfs(scantab, scannos=[], smooth=1, tsysval=0.0, tauval=0.0, tcalval=0.0, verify=False): 00622 """ 00623 Calibrate GBT frequency switched data. 00624 Adopted from GBTIDL getfs. 00625 Currently calfs identify the scans as frequency switched data if source 00626 type enum is fson and fsoff. The data must contains 'CAL' signal 00627 on/off in each integration. To identify 'CAL' on state, the source type 00628 enum of foncal and foffcal need to be present. 00629 00630 Parameters: 00631 scantab: scantable 00632 scannos: list of scan numbers 00633 smooth: optional box smoothing order for the reference 00634 (default is 1 = no smoothing) 00635 tsysval: optional user specified Tsys (default is 0.0, 00636 use Tsys in the data) 00637 tauval: optional user specified Tau 00638 verify: Verify calibration if true 00639 """ 00640 varlist = vars() 00641 from asap._asap import stmath 00642 from asap._asap import srctype 00643 stm = stmath() 00644 stm._setinsitu(False) 00645 00646 # check = scantab.get_scan('*_fs*') 00647 # if check is None: 00648 # msg = "The input data appear to contain no Nod observing mode data." 00649 # raise TypeError(msg) 00650 s = scantab.get_scan(scannos) 00651 del scantab 00652 00653 resspec = scantable(stm._dofs(s, scannos, smooth, tsysval,tauval,tcalval)) 00654 ### 00655 if verify: 00656 # get data 00657 ssub = s.get_scan(scannos) 00658 #ssubon = ssub.get_scan('*calon') 00659 #ssuboff = ssub.get_scan('*[^calon]') 00660 sel = selector() 00661 sel.set_types( [srctype.foncal,srctype.foffcal] ) 00662 ssub.set_selection( sel ) 00663 ssubon = ssub.copy() 00664 ssub.set_selection() 00665 sel.reset() 00666 sel.set_types( [srctype.fson,srctype.fsoff] ) 00667 ssub.set_selection( sel ) 00668 ssuboff = ssub.copy() 00669 ssub.set_selection() 00670 sel.reset() 00671 import numpy 00672 precal={} 00673 postcal=[] 00674 keys=['fs','fs_calon','fsr','fsr_calon'] 00675 types=[srctype.fson,srctype.foncal,srctype.fsoff,srctype.foffcal] 00676 ifnos=list(ssub.getifnos()) 00677 polnos=list(ssub.getpolnos()) 00678 for i in range(2): 00679 #ss=ssuboff.get_scan('*'+keys[2*i]) 00680 ll=[] 00681 for j in range(len(ifnos)): 00682 for k in range(len(polnos)): 00683 sel.set_ifs(ifnos[j]) 00684 sel.set_polarizations(polnos[k]) 00685 sel.set_types(types[2*i]) 00686 try: 00687 #ss.set_selection(sel) 00688 ssuboff.set_selection(sel) 00689 except: 00690 continue 00691 ll.append(numpy.array(ss._getspectrum(0))) 00692 sel.reset() 00693 #ss.set_selection() 00694 ssuboff.set_selection() 00695 precal[keys[2*i]]=ll 00696 #del ss 00697 #ss=ssubon.get_scan('*'+keys[2*i+1]) 00698 ll=[] 00699 for j in range(len(ifnos)): 00700 for k in range(len(polnos)): 00701 sel.set_ifs(ifnos[j]) 00702 sel.set_polarizations(polnos[k]) 00703 sel.set_types(types[2*i+1]) 00704 try: 00705 #ss.set_selection(sel) 00706 ssubon.set_selection(sel) 00707 except: 00708 continue 00709 ll.append(numpy.array(ss._getspectrum(0))) 00710 sel.reset() 00711 #ss.set_selection() 00712 ssubon.set_selection() 00713 precal[keys[2*i+1]]=ll 00714 #del ss 00715 #sig=resspec.get_scan('*_fs') 00716 #ref=resspec.get_scan('*_fsr') 00717 sel.set_types( srctype.fson ) 00718 resspec.set_selection( sel ) 00719 sig=resspec.copy() 00720 resspec.set_selection() 00721 sel.reset() 00722 sel.set_type( srctype.fsoff ) 00723 resspec.set_selection( sel ) 00724 ref=resspec.copy() 00725 resspec.set_selection() 00726 sel.reset() 00727 for k in range(len(polnos)): 00728 for j in range(len(ifnos)): 00729 sel.set_ifs(ifnos[j]) 00730 sel.set_polarizations(polnos[k]) 00731 try: 00732 sig.set_selection(sel) 00733 postcal.append(numpy.array(sig._getspectrum(0))) 00734 except: 00735 ref.set_selection(sel) 00736 postcal.append(numpy.array(ref._getspectrum(0))) 00737 sel.reset() 00738 resspec.set_selection() 00739 del sel 00740 # plot 00741 asaplog.post() 00742 asaplog.push('Plot only first spectrum for each [if,pol] pairs to verify calibration.') 00743 asaplog.post('WARN') 00744 p=new_asaplot() 00745 rcp('lines', linewidth=1) 00746 #nr=min(6,len(ifnos)*len(polnos)) 00747 nr=len(ifnos)/2*len(polnos) 00748 titles=[] 00749 btics=[] 00750 if nr>3: 00751 asaplog.post() 00752 asaplog.push('Only first 3 [if,pol] pairs are plotted.') 00753 asaplog.post('WARN') 00754 nr=3 00755 p.set_panels(rows=nr,cols=3,nplots=3*nr,ganged=False) 00756 for i in range(3*nr): 00757 b=False 00758 if i >= 3*nr-3: 00759 b=True 00760 btics.append(b) 00761 for i in range(nr): 00762 p.subplot(3*i) 00763 p.color=0 00764 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)]) 00765 titles.append(title) 00766 #p.set_axes('title',title,fontsize=40) 00767 ymin=1.0e100 00768 ymax=-1.0e100 00769 nchan=s.nchan(ifnos[2*int(i/len(polnos))]) 00770 edge=int(nchan*0.01) 00771 for j in range(4): 00772 spmin=min(precal[keys[j]][i][edge:nchan-edge]) 00773 spmax=max(precal[keys[j]][i][edge:nchan-edge]) 00774 ymin=min(ymin,spmin) 00775 ymax=max(ymax,spmax) 00776 for j in range(4): 00777 if i==0: 00778 p.set_line(label=keys[j]) 00779 else: 00780 p.legend() 00781 p.plot(precal[keys[j]][i]) 00782 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax)) 00783 if not btics[3*i]: 00784 p.axes.set_xticks([]) 00785 p.subplot(3*i+1) 00786 p.color=0 00787 title='sig data IF%s POL%s' % (ifnos[2*int(i/len(polnos))],polnos[i%len(polnos)]) 00788 titles.append(title) 00789 #p.set_axes('title',title) 00790 p.legend() 00791 ymin=postcal[2*i][edge:nchan-edge].min() 00792 ymax=postcal[2*i][edge:nchan-edge].max() 00793 p.plot(postcal[2*i]) 00794 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax)) 00795 if not btics[3*i+1]: 00796 p.axes.set_xticks([]) 00797 p.subplot(3*i+2) 00798 p.color=0 00799 title='ref data IF%s POL%s' % (ifnos[2*int(i/len(polnos))+1],polnos[i%len(polnos)]) 00800 titles.append(title) 00801 #p.set_axes('title',title) 00802 p.legend() 00803 ymin=postcal[2*i+1][edge:nchan-edge].min() 00804 ymax=postcal[2*i+1][edge:nchan-edge].max() 00805 p.plot(postcal[2*i+1]) 00806 p.axes.set_ylim(ymin-0.1*abs(ymin),ymax+0.1*abs(ymax)) 00807 if not btics[3*i+2]: 00808 p.axes.set_xticks([]) 00809 for i in range(3*nr): 00810 p.subplot(i) 00811 p.set_axes('title',titles[i],fontsize='medium') 00812 x=raw_input('Accept calibration ([y]/n): ' ) 00813 if x.upper() == 'N': 00814 p.quit() 00815 del p 00816 return scabtab 00817 p.quit() 00818 del p 00819 ### 00820 resspec._add_history("calfs",varlist) 00821 return resspec 00822 00823 @asaplog_post_dec 00824 def merge(*args): 00825 """ 00826 Merge a list of scanatables, or comma-sperated scantables into one 00827 scnatble. 00828 Parameters: 00829 A list [scan1, scan2] or scan1, scan2. 00830 Example: 00831 myscans = [scan1, scan2] 00832 allscans = merge(myscans) 00833 # or equivalent 00834 sameallscans = merge(scan1, scan2) 00835 """ 00836 varlist = vars() 00837 if isinstance(args[0],list): 00838 lst = tuple(args[0]) 00839 elif isinstance(args[0],tuple): 00840 lst = args[0] 00841 else: 00842 lst = tuple(args) 00843 varlist["args"] = "%d scantables" % len(lst) 00844 # need special formatting her for history... 00845 from asap._asap import stmath 00846 stm = stmath() 00847 for s in lst: 00848 if not isinstance(s,scantable): 00849 msg = "Please give a list of scantables" 00850 raise TypeError(msg) 00851 s = scantable(stm._merge(lst)) 00852 s._add_history("merge", varlist) 00853 return s 00854 00855 @asaplog_post_dec 00856 def calibrate( scantab, scannos=[], calmode='none', verify=None ): 00857 """ 00858 Calibrate data. 00859 00860 Parameters: 00861 scantab: scantable 00862 scannos: list of scan number 00863 calmode: calibration mode 00864 verify: verify calibration 00865 """ 00866 import re 00867 antname = scantab.get_antennaname() 00868 if ( calmode == 'nod' ): 00869 asaplog.push( 'Calibrating nod data.' ) 00870 scal = calnod( scantab, scannos=scannos, verify=verify ) 00871 elif ( calmode == 'quotient' ): 00872 asaplog.push( 'Calibrating using quotient.' ) 00873 scal = scantab.auto_quotient( verify=verify ) 00874 elif ( calmode == 'ps' ): 00875 asaplog.push( 'Calibrating %s position-switched data.' % antname ) 00876 if ( antname.find( 'APEX' ) != -1 ): 00877 scal = apexcal( scantab, scannos, calmode, verify ) 00878 elif ( antname.find( 'ALMA' ) != -1 or antname.find( 'OSF' ) != -1 00879 or re.match('DV[0-9][0-9]$',antname) is not None 00880 or re.match('PM[0-9][0-9]$',antname) is not None 00881 or re.match('CM[0-9][0-9]$',antname) is not None 00882 or re.match('DA[0-9][0-9]$',antname) is not None ): 00883 scal = almacal( scantab, scannos, calmode, verify ) 00884 else: 00885 scal = calps( scantab, scannos=scannos, verify=verify ) 00886 elif ( calmode == 'fs' or calmode == 'fsotf' ): 00887 asaplog.push( 'Calibrating %s frequency-switched data.' % antname ) 00888 if ( antname.find( 'APEX' ) != -1 ): 00889 scal = apexcal( scantab, scannos, calmode, verify ) 00890 elif ( antname.find( 'ALMA' ) != -1 or antname.find( 'OSF' ) != -1 ): 00891 scal = almacal( scantab, scannos, calmode, verify ) 00892 else: 00893 scal = calfs( scantab, scannos=scannos, verify=verify ) 00894 elif ( calmode == 'otf' ): 00895 asaplog.push( 'Calibrating %s On-The-Fly data.' % antname ) 00896 scal = almacal( scantab, scannos, calmode, verify ) 00897 else: 00898 asaplog.push( 'No calibration.' ) 00899 scal = scantab.copy() 00900 00901 return scal 00902 00903 def apexcal( scantab, scannos=[], calmode='none', verify=False ): 00904 """ 00905 Calibrate APEX data 00906 00907 Parameters: 00908 scantab: scantable 00909 scannos: list of scan number 00910 calmode: calibration mode 00911 00912 verify: verify calibration 00913 """ 00914 from asap._asap import stmath 00915 stm = stmath() 00916 antname = scantab.get_antennaname() 00917 selection=selector() 00918 selection.set_scans(scannos) 00919 orig = scantab.get_selection() 00920 scantab.set_selection(orig+selection) 00921 ## ssub = scantab.get_scan( scannos ) 00922 ## scal = scantable( stm.cwcal( ssub, calmode, antname ) ) 00923 scal = scantable( stm.cwcal( scantab, calmode, antname ) ) 00924 scantab.set_selection(orig) 00925 return scal 00926 00927 def almacal( scantab, scannos=[], calmode='none', verify=False ): 00928 """ 00929 Calibrate ALMA data 00930 00931 Parameters: 00932 scantab: scantable 00933 scannos: list of scan number 00934 calmode: calibration mode 00935 00936 verify: verify calibration 00937 """ 00938 from asap._asap import stmath 00939 stm = stmath() 00940 selection=selector() 00941 selection.set_scans(scannos) 00942 orig = scantab.get_selection() 00943 scantab.set_selection(orig+selection) 00944 ## ssub = scantab.get_scan( scannos ) 00945 ## scal = scantable( stm.almacal( ssub, calmode ) ) 00946 scal = scantable( stm.almacal( scantab, calmode ) ) 00947 scantab.set_selection(orig) 00948 return scal 00949 00950 @asaplog_post_dec 00951 def splitant(filename, outprefix='',overwrite=False): 00952 """ 00953 Split Measurement set by antenna name, save data as a scantables, 00954 and return a list of filename. 00955 Notice this method can only be available from CASA. 00956 Prameter 00957 filename: the name of Measurement set to be read. 00958 outprefix: the prefix of output scantable name. 00959 the names of output scantable will be 00960 outprefix.antenna1, outprefix.antenna2, .... 00961 If not specified, outprefix = filename is assumed. 00962 overwrite If the file should be overwritten if it exists. 00963 The default False is to return with warning 00964 without writing the output. USE WITH CARE. 00965 00966 """ 00967 # Import the table toolkit from CASA 00968 from casac import casac 00969 from asap.scantable import is_ms 00970 tb = casac.table() 00971 # Check the input filename 00972 if isinstance(filename, str): 00973 import os.path 00974 filename = os.path.expandvars(filename) 00975 filename = os.path.expanduser(filename) 00976 if not os.path.exists(filename): 00977 s = "File '%s' not found." % (filename) 00978 raise IOError(s) 00979 # check if input file is MS 00980 #if not os.path.isdir(filename) \ 00981 # or not os.path.exists(filename+'/ANTENNA') \ 00982 # or not os.path.exists(filename+'/table.f1'): 00983 if not is_ms(filename): 00984 s = "File '%s' is not a Measurement set." % (filename) 00985 raise IOError(s) 00986 else: 00987 s = "The filename should be string. " 00988 raise TypeError(s) 00989 # Check out put file name 00990 outname='' 00991 if len(outprefix) > 0: prefix=outprefix+'.' 00992 else: 00993 prefix=filename.rstrip('/') 00994 # Now do the actual splitting. 00995 outfiles=[] 00996 tb.open(tablename=filename,nomodify=True) 00997 ant1=tb.getcol('ANTENNA1',0,-1,1) 00998 #anttab=tb.getkeyword('ANTENNA').split()[-1] 00999 anttab=tb.getkeyword('ANTENNA').lstrip('Table: ') 01000 tb.close() 01001 #tb.open(tablename=filename+'/ANTENNA',nomodify=True) 01002 tb.open(tablename=anttab,nomodify=True) 01003 nant=tb.nrows() 01004 antnames=tb.getcol('NAME',0,nant,1) 01005 tb.close() 01006 tmpname='asapmath.splitant.tmp' 01007 for antid in set(ant1): 01008 tb.open(tablename=filename,nomodify=True) 01009 tbsel=tb.query('ANTENNA1 == %s && ANTENNA2 == %s'%(antid,antid),tmpname) 01010 scan=scantable(tmpname,average=False,getpt=True,antenna=int(antid)) 01011 outname=prefix+antnames[antid]+'.asap' 01012 scan.save(outname,format='ASAP',overwrite=overwrite) 01013 tbsel.close() 01014 tb.close() 01015 del tbsel 01016 del scan 01017 outfiles.append(outname) 01018 os.system('rm -rf '+tmpname) 01019 del tb 01020 return outfiles 01021 01022 @asaplog_post_dec 01023 def _array2dOp( scan, value, mode="ADD", tsys=False, insitu=None): 01024 """ 01025 This function is workaround on the basic operation of scantable 01026 with 2 dimensional float list. 01027 01028 scan: scantable operand 01029 value: float list operand 01030 mode: operation mode (ADD, SUB, MUL, DIV) 01031 tsys: if True, operate tsys as well 01032 insitu: if False, a new scantable is returned. 01033 Otherwise, the array operation is done in-sitsu. 01034 """ 01035 if insitu is None: insitu = rcParams['insitu'] 01036 nrow = scan.nrow() 01037 s = None 01038 from asap._asap import stmath 01039 stm = stmath() 01040 stm._setinsitu(insitu) 01041 if len( value ) == 1: 01042 s = scantable( stm._arrayop( scan, value[0], mode, tsys ) ) 01043 elif len( value ) != nrow: 01044 raise ValueError( 'len(value) must be 1 or conform to scan.nrow()' ) 01045 else: 01046 from asap._asap import stmath 01047 if not insitu: 01048 s = scan.copy() 01049 else: 01050 s = scan 01051 # insitu must be True as we go row by row on the same data 01052 stm._setinsitu( True ) 01053 basesel = s.get_selection() 01054 # generate a new selector object based on basesel 01055 sel = selector(basesel) 01056 for irow in range( nrow ): 01057 sel.set_rows( irow ) 01058 s.set_selection( sel ) 01059 if len( value[irow] ) == 1: 01060 stm._unaryop( s, value[irow][0], mode, tsys ) 01061 else: 01062 #stm._arrayop( s, value[irow], mode, tsys, 'channel' ) 01063 stm._arrayop( s, value[irow], mode, tsys ) 01064 s.set_selection(basesel) 01065 return s
1.8.0
