Changeset 1789 for branches/mergetest/python/scantable.py
- Timestamp:
- 07/30/10 16:59:56 (14 years ago)
- Location:
- branches/mergetest/python
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
branches/mergetest/python
-
Property
svn:mergeinfo
set to
/branches/alma/python merged eligible
-
Property
svn:mergeinfo
set to
-
branches/mergetest/python/scantable.py
r1731 r1789 35 35 The ASAP container for scans 36 36 """ 37 @print_log_dec 38 def __init__(self, filename, average=None, unit=None, parallactify=None): 37 38 @print_log_dec 39 def __init__(self, filename, average=None, unit=None, getpt=None, antenna=None, parallactify=None): 39 40 """ 40 41 Create a scantable from a saved one or make a reference … … 54 55 (input rpfits/sdfits/ms) or replaces the value 55 56 in existing scantables 57 getpt: for MeasurementSet input data only: 58 If True, all pointing data are filled. 59 The deafult is False, which makes time to load 60 the MS data faster in some cases. 61 antenna: Antenna selection. integer (id) or string (name 62 or id). 56 63 parallactify: Indcicate that the data had been parallatified. 57 64 Default is taken form rc file. … … 59 66 if average is None: 60 67 average = rcParams['scantable.autoaverage'] 68 if getpt is None: 69 getpt = True 70 if antenna is None: 71 antenna = '' 72 elif type(antenna) == int: 73 antenna = '%s'%antenna 74 elif type(antenna) == list: 75 tmpstr = '' 76 for i in range( len(antenna) ): 77 if type(antenna[i]) == int: 78 tmpstr = tmpstr + ('%s,'%(antenna[i])) 79 elif type(antenna[i]) == str: 80 tmpstr=tmpstr+antenna[i]+',' 81 else: 82 asaplog.push('Bad antenna selection.') 83 print_log('ERROR') 84 return 85 antenna = tmpstr.rstrip(',') 61 86 parallactify = parallactify or rcParams['scantable.parallactify'] 62 87 varlist = vars() 63 88 from asap._asap import stmath 64 self._math = stmath( )89 self._math = stmath( rcParams['insitu'] ) 65 90 if isinstance(filename, Scantable): 66 91 Scantable.__init__(self, filename) … … 73 98 if rcParams['verbose']: 74 99 asaplog.push(s) 100 print_log('ERROR') 75 101 return 76 102 raise IOError(s) … … 80 106 if unit is not None: 81 107 self.set_fluxunit(unit) 82 self.set_freqframe(rcParams['scantable.freqframe']) 108 # do not reset to the default freqframe 109 #self.set_freqframe(rcParams['scantable.freqframe']) 110 elif os.path.isdir(filename) \ 111 and not os.path.exists(filename+'/table.f1'): 112 msg = "The given file '%s'is not a valid " \ 113 "asap table." % (filename) 114 if rcParams['verbose']: 115 #print msg 116 asaplog.push( msg ) 117 print_log( 'ERROR' ) 118 return 119 else: 120 raise IOError(msg) 83 121 else: 84 self._fill([filename], unit, average )122 self._fill([filename], unit, average, getpt, antenna) 85 123 elif (isinstance(filename, list) or isinstance(filename, tuple)) \ 86 124 and isinstance(filename[-1], str): 87 self._fill(filename, unit, average )125 self._fill(filename, unit, average, getpt, antenna) 88 126 self.parallactify(parallactify) 89 127 self._add_history("scantable", varlist) 128 print_log() 90 129 91 130 @print_log_dec … … 126 165 msg = "File %s exists." % name 127 166 if rcParams['verbose']: 128 print msg 167 #print msg 168 asaplog.push( msg ) 169 print_log( 'ERROR' ) 129 170 return 130 171 else: … … 137 178 writer = stw(format2) 138 179 writer.write(self, name) 180 print_log() 139 181 return 140 182 … … 164 206 if not _is_valid(scanid): 165 207 if rcParams['verbose']: 166 print "Please specify a scanno to drop from the scantable" 208 #print "Please specify a scanno to drop from the scantable" 209 asaplog.push( 'Please specify a scanno to drop from the scantable' ) 210 print_log( 'ERROR' ) 167 211 return 168 212 else: … … 176 220 except ValueError: 177 221 if rcParams['verbose']: 178 print "Couldn't find any match." 222 #print "Couldn't find any match." 223 print_log() 224 asaplog.push( "Couldn't find any match." ) 225 print_log( 'ERROR' ) 179 226 return 180 227 else: raise … … 184 231 except RuntimeError: 185 232 if rcParams['verbose']: 186 print "Couldn't find any match." 233 #print "Couldn't find any match." 234 print_log() 235 asaplog.push( "Couldn't find any match." ) 236 print_log( 'ERROR' ) 187 237 else: 188 238 raise … … 215 265 if scanid is None: 216 266 if rcParams['verbose']: 217 print "Please specify a scan no or name to " \ 218 "retrieve from the scantable" 267 #print "Please specify a scan no or name to " \ 268 # "retrieve from the scantable" 269 asaplog.push( 'Please specify a scan no or name to retrieve from the scantable' ) 270 print_log( 'ERROR' ) 219 271 return 220 272 else: … … 236 288 msg = "Illegal scanid type, use 'int' or 'list' if ints." 237 289 if rcParams['verbose']: 238 print msg 290 #print msg 291 asaplog.push( msg ) 292 print_log( 'ERROR' ) 239 293 else: 240 294 raise TypeError(msg) 241 295 except RuntimeError: 242 if rcParams['verbose']: print "Couldn't find any match." 296 if rcParams['verbose']: 297 #print "Couldn't find any match." 298 print_log() 299 asaplog.push( "Couldn't find any match." ) 300 print_log( 'ERROR' ) 243 301 else: raise 244 302 … … 266 324 msg = "Illegal file name '%s'." % (filename) 267 325 if rcParams['verbose']: 268 print msg 326 #print msg 327 asaplog.push( msg ) 328 print_log( 'ERROR' ) 269 329 else: 270 330 raise IOError(msg) … … 367 427 self._setselection(selection) 368 428 369 def stats(self, stat='stddev', mask=None): 429 def get_row(self, row=0, insitu=None): 430 """ 431 Select a row in the scantable. 432 Return a scantable with single row. 433 Parameters: 434 row: row no of integration, default is 0. 435 insitu: if False a new scantable is returned. 436 Otherwise, the scaling is done in-situ 437 The default is taken from .asaprc (False) 438 """ 439 if insitu is None: insitu = rcParams['insitu'] 440 if not insitu: 441 workscan = self.copy() 442 else: 443 workscan = self 444 # Select a row 445 sel=selector() 446 sel.set_scans([workscan.getscan(row)]) 447 sel.set_cycles([workscan.getcycle(row)]) 448 sel.set_beams([workscan.getbeam(row)]) 449 sel.set_ifs([workscan.getif(row)]) 450 sel.set_polarisations([workscan.getpol(row)]) 451 sel.set_name(workscan._getsourcename(row)) 452 workscan.set_selection(sel) 453 if not workscan.nrow() == 1: 454 msg = "Cloud not identify single row. %d rows selected."%(workscan.nrow()) 455 raise RuntimeError(msg) 456 del sel 457 if insitu: 458 self._assign(workscan) 459 else: 460 return workscan 461 462 #def stats(self, stat='stddev', mask=None): 463 def stats(self, stat='stddev', mask=None, form='3.3f'): 370 464 """ 371 465 Determine the specified statistic of the current beam/if/pol … … 373 467 channels should be excluded. 374 468 Parameters: 375 stat: 'min', 'max', ' sumsq', 'sum', 'mean'376 ' var', 'stddev', 'avdev', 'rms', 'median'469 stat: 'min', 'max', 'min_abc', 'max_abc', 'sumsq', 'sum', 470 'mean', 'var', 'stddev', 'avdev', 'rms', 'median' 377 471 mask: an optional mask specifying where the statistic 378 472 should be determined. 473 form: format string to print statistic values 379 474 Example: 380 475 scan.set_unit('channel') … … 387 482 "number of channels. Please use setselection() " 388 483 "to select individual IFs") 389 390 statvals = self._math._stats(self, mask, stat) 391 def cb(i): 392 return statvals[i] 393 394 return self._row_callback(cb, stat) 484 rtnabc = False 485 if stat.lower().endswith('_abc'): rtnabc = True 486 getchan = False 487 if stat.lower().startswith('min') or stat.lower().startswith('max'): 488 chan = self._math._minmaxchan(self, mask, stat) 489 getchan = True 490 statvals = [] 491 if not rtnabc: statvals = self._math._stats(self, mask, stat) 492 493 #def cb(i): 494 # return statvals[i] 495 496 #return self._row_callback(cb, stat) 497 498 label=stat 499 #callback=cb 500 out = "" 501 #outvec = [] 502 sep = '-'*50 503 for i in range(self.nrow()): 504 refstr = '' 505 statunit= '' 506 if getchan: 507 qx, qy = self.chan2data(rowno=i, chan=chan[i]) 508 if rtnabc: 509 statvals.append(qx['value']) 510 refstr = ('(value: %'+form) % (qy['value'])+' ['+qy['unit']+'])' 511 statunit= '['+qx['unit']+']' 512 else: 513 refstr = ('(@ %'+form) % (qx['value'])+' ['+qx['unit']+'])' 514 515 tm = self._gettime(i) 516 src = self._getsourcename(i) 517 out += 'Scan[%d] (%s) ' % (self.getscan(i), src) 518 out += 'Time[%s]:\n' % (tm) 519 if self.nbeam(-1) > 1: 520 out += ' Beam[%d] ' % (self.getbeam(i)) 521 if self.nif(-1) > 1: out += ' IF[%d] ' % (self.getif(i)) 522 if self.npol(-1) > 1: out += ' Pol[%d] ' % (self.getpol(i)) 523 #outvec.append(callback(i)) 524 #out += ('= %'+form) % (outvec[i]) +' '+refstr+'\n' 525 out += ('= %'+form) % (statvals[i]) +' '+refstr+'\n' 526 out += sep+"\n" 527 528 if rcParams['verbose']: 529 import os 530 if os.environ.has_key( 'USER' ): 531 usr=os.environ['USER'] 532 else: 533 import commands 534 usr=commands.getoutput( 'whoami' ) 535 tmpfile='/tmp/tmp_'+usr+'_casapy_asap_scantable_stats' 536 f=open(tmpfile,'w') 537 print >> f, sep 538 print >> f, ' %s %s' % (label, statunit) 539 print >> f, sep 540 print >> f, out 541 f.close() 542 f=open(tmpfile,'r') 543 x=f.readlines() 544 f.close() 545 blanc='' 546 asaplog.push(blanc.join(x), False) 547 #for xx in x: 548 # asaplog.push( xx, False ) 549 print_log() 550 return statvals 551 552 def chan2data(self, rowno=0, chan=0): 553 """ 554 Returns channel/frequency/velocity and spectral value 555 at an arbitrary row and channel in the scantable. 556 Parameters: 557 rowno: a row number in the scantable. Default is the 558 first row, i.e. rowno=0 559 chan: a channel in the scantable. Default is the first 560 channel, i.e. pos=0 561 """ 562 if isinstance(rowno, int) and isinstance(chan, int): 563 qx = {'unit': self.get_unit(), 564 'value': self._getabcissa(rowno)[chan]} 565 qy = {'unit': self.get_fluxunit(), 566 'value': self._getspectrum(rowno)[chan]} 567 return qx, qy 395 568 396 569 def stddev(self, mask=None): … … 462 635 out += sep+'\n' 463 636 if rcParams['verbose']: 464 print sep 465 print " %s" % (label) 466 print sep 467 print out 637 asaplog.push(sep) 638 asaplog.push(" %s" % (label)) 639 asaplog.push(sep) 640 asaplog.push(out) 641 print_log() 468 642 return outvec 469 643 … … 474 648 return [callback(i) for i in range(self.nrow())] 475 649 else: 476 if 0 <= row < self.nrow():650 if 0 <= row < self.nrow(): 477 651 return callback(row) 478 652 … … 604 778 self._setInstrument(instr) 605 779 self._add_history("set_instument", vars()) 780 print_log() 606 781 607 782 @print_log_dec … … 614 789 self._setfeedtype(feedtype) 615 790 self._add_history("set_feedtype", vars()) 791 print_log() 616 792 617 793 @print_log_dec … … 627 803 self._setcoordinfo(inf) 628 804 self._add_history("set_doppler", vars()) 805 print_log() 629 806 630 807 @print_log_dec … … 634 811 Parameters: 635 812 frame: an optional frame type, default 'LSRK'. Valid frames are: 636 ' REST', 'TOPO', 'LSRD', 'LSRK', 'BARY',813 'TOPO', 'LSRD', 'LSRK', 'BARY', 637 814 'GEO', 'GALACTO', 'LGROUP', 'CMB' 638 815 Examples: … … 641 818 frame = frame or rcParams['scantable.freqframe'] 642 819 varlist = vars() 643 valid = ['REST', 'TOPO', 'LSRD', 'LSRK', 'BARY', \ 820 # "REST" is not implemented in casacore 821 #valid = ['REST', 'TOPO', 'LSRD', 'LSRK', 'BARY', \ 822 # 'GEO', 'GALACTO', 'LGROUP', 'CMB'] 823 valid = ['TOPO', 'LSRD', 'LSRK', 'BARY', \ 644 824 'GEO', 'GALACTO', 'LGROUP', 'CMB'] 645 825 … … 652 832 msg = "Please specify a valid freq type. Valid types are:\n", valid 653 833 if rcParams['verbose']: 654 print msg 834 #print msg 835 asaplog.push( msg ) 836 print_log( 'ERROR' ) 655 837 else: 656 838 raise TypeError(msg) 839 print_log() 657 840 658 841 def set_dirframe(self, frame=""): … … 670 853 except RuntimeError, msg: 671 854 if rcParams['verbose']: 672 print msg 855 #print msg 856 print_log() 857 asaplog.push( str(msg) ) 858 print_log( 'ERROR' ) 673 859 else: 674 860 raise … … 698 884 abc = self._getabcissa(rowno) 699 885 lbl = self._getabcissalabel(rowno) 886 print_log() 700 887 return abc, lbl 701 888 702 def flag(self, mask=None ):889 def flag(self, mask=None, unflag=False): 703 890 """ 704 891 Flag the selected data using an optional channel mask. … … 706 893 mask: an optional channel mask, created with create_mask. Default 707 894 (no mask) is all channels. 895 unflag: if True, unflag the data 708 896 """ 709 897 varlist = vars() 710 898 mask = mask or [] 711 899 try: 712 self._flag(mask )900 self._flag(mask, unflag) 713 901 except RuntimeError, msg: 714 902 if rcParams['verbose']: 715 print msg 903 #print msg 904 print_log() 905 asaplog.push( str(msg) ) 906 print_log( 'ERROR' ) 716 907 return 717 908 else: raise 718 909 self._add_history("flag", varlist) 719 910 911 def flag_row(self, rows=[], unflag=False): 912 """ 913 Flag the selected data in row-based manner. 914 Parameters: 915 rows: list of row numbers to be flagged. Default is no row (must be explicitly specified to execute row-based flagging). 916 unflag: if True, unflag the data. 917 """ 918 varlist = vars() 919 try: 920 self._flag_row(rows, unflag) 921 except RuntimeError, msg: 922 if rcParams['verbose']: 923 print_log() 924 asaplog.push( str(msg) ) 925 print_log('ERROR') 926 return 927 else: raise 928 self._add_history("flag_row", varlist) 929 930 def clip(self, uthres=None, dthres=None, clipoutside=True, unflag=False): 931 """ 932 Flag the selected data outside a specified range (in channel-base) 933 Parameters: 934 uthres: upper threshold. 935 dthres: lower threshold 936 clipoutside: True for flagging data outside the range [dthres:uthres]. 937 False for glagging data inside the range. 938 unflag : if True, unflag the data. 939 """ 940 varlist = vars() 941 try: 942 self._clip(uthres, dthres, clipoutside, unflag) 943 except RuntimeError, msg: 944 if rcParams['verbose']: 945 print_log() 946 asaplog.push(str(msg)) 947 print_log('ERROR') 948 return 949 else: raise 950 self._add_history("clip", varlist) 951 720 952 @print_log_dec 721 953 def lag_flag(self, start, end, unit="MHz", insitu=None): 954 #def lag_flag(self, frequency, width=0.0, unit="GHz", insitu=None): 722 955 """ 723 956 Flag the data in 'lag' space by providing a frequency to remove. … … 748 981 except RuntimeError, msg: 749 982 if rcParams['verbose']: 750 print msg 983 #print msg 984 print_log() 985 asaplog.push( str(msg) ) 986 print_log( 'ERROR' ) 751 987 return 752 988 else: raise 753 989 s._add_history("lag_flag", varlist) 990 print_log() 754 991 if insitu: 755 992 self._assign(s) … … 819 1056 if kwargs.get('invert'): 820 1057 msk = mask_not(msk) 1058 print_log() 821 1059 return msk 822 1060 823 def get_restfreqs(self): 1061 def get_masklist(self, mask=None, row=0): 1062 """ 1063 Compute and return a list of mask windows, [min, max]. 1064 Parameters: 1065 mask: channel mask, created with create_mask. 1066 row: calcutate the masklist using the specified row 1067 for unit conversions, default is row=0 1068 only necessary if frequency varies over rows. 1069 Returns: 1070 [min, max], [min2, max2], ... 1071 Pairs of start/end points (inclusive)specifying 1072 the masked regions 1073 """ 1074 if not (isinstance(mask,list) or isinstance(mask, tuple)): 1075 raise TypeError("The mask should be list or tuple.") 1076 if len(mask) < 2: 1077 raise TypeError("The mask elements should be > 1") 1078 if self.nchan() != len(mask): 1079 msg = "Number of channels in scantable != number of mask elements" 1080 raise TypeError(msg) 1081 data = self._getabcissa(row) 1082 u = self._getcoordinfo()[0] 1083 if rcParams['verbose']: 1084 if u == "": u = "channel" 1085 msg = "The current mask window unit is %s" % u 1086 i = self._check_ifs() 1087 if not i: 1088 msg += "\nThis mask is only valid for IF=%d" % (self.getif(i)) 1089 asaplog.push(msg) 1090 masklist=[] 1091 ist, ien = None, None 1092 ist, ien=self.get_mask_indices(mask) 1093 if ist is not None and ien is not None: 1094 for i in xrange(len(ist)): 1095 range=[data[ist[i]],data[ien[i]]] 1096 range.sort() 1097 masklist.append([range[0],range[1]]) 1098 return masklist 1099 1100 def get_mask_indices(self, mask=None): 1101 """ 1102 Compute and Return lists of mask start indices and mask end indices. 1103 Parameters: 1104 mask: channel mask, created with create_mask. 1105 Returns: 1106 List of mask start indices and that of mask end indices, 1107 i.e., [istart1,istart2,....], [iend1,iend2,....]. 1108 """ 1109 if not (isinstance(mask,list) or isinstance(mask, tuple)): 1110 raise TypeError("The mask should be list or tuple.") 1111 if len(mask) < 2: 1112 raise TypeError("The mask elements should be > 1") 1113 istart=[] 1114 iend=[] 1115 if mask[0]: istart.append(0) 1116 for i in range(len(mask)-1): 1117 if not mask[i] and mask[i+1]: 1118 istart.append(i+1) 1119 elif mask[i] and not mask[i+1]: 1120 iend.append(i) 1121 if mask[len(mask)-1]: iend.append(len(mask)-1) 1122 if len(istart) != len(iend): 1123 raise RuntimeError("Numbers of mask start != mask end.") 1124 for i in range(len(istart)): 1125 if istart[i] > iend[i]: 1126 raise RuntimeError("Mask start index > mask end index") 1127 break 1128 return istart,iend 1129 1130 # def get_restfreqs(self): 1131 # """ 1132 # Get the restfrequency(s) stored in this scantable. 1133 # The return value(s) are always of unit 'Hz' 1134 # Parameters: 1135 # none 1136 # Returns: 1137 # a list of doubles 1138 # """ 1139 # return list(self._getrestfreqs()) 1140 1141 def get_restfreqs(self, ids=None): 824 1142 """ 825 1143 Get the restfrequency(s) stored in this scantable. 826 1144 The return value(s) are always of unit 'Hz' 827 1145 Parameters: 828 none 1146 ids: (optional) a list of MOLECULE_ID for that restfrequency(s) to 1147 be retrieved 829 1148 Returns: 830 a list of doubles 831 """ 832 return list(self._getrestfreqs()) 833 1149 dictionary containing ids and a list of doubles for each id 1150 """ 1151 if ids is None: 1152 rfreqs={} 1153 idlist = self.getmolnos() 1154 for i in idlist: 1155 rfreqs[i]=list(self._getrestfreqs(i)) 1156 return rfreqs 1157 else: 1158 if type(ids)==list or type(ids)==tuple: 1159 rfreqs={} 1160 for i in ids: 1161 rfreqs[i]=list(self._getrestfreqs(i)) 1162 return rfreqs 1163 else: 1164 return list(self._getrestfreqs(ids)) 1165 #return list(self._getrestfreqs(ids)) 834 1166 835 1167 def set_restfreqs(self, freqs=None, unit='Hz'): 836 1168 """ 1169 ********NEED TO BE UPDATED begin************ 837 1170 Set or replace the restfrequency specified and 838 1171 If the 'freqs' argument holds a scalar, … … 846 1179 E.g. 'freqs=[1e9, 2e9]' would mean IF 0 gets restfreq 1e9 and 847 1180 IF 1 gets restfreq 2e9. 1181 ********NEED TO BE UPDATED end************ 848 1182 You can also specify the frequencies via a linecatalog. 849 1183 … … 853 1187 854 1188 Example: 855 # set the given restfrequency for the whole table1189 # set the given restfrequency for the all currently selected IFs 856 1190 scan.set_restfreqs(freqs=1.4e9) 857 # If thee number of IFs in the data is >= 2 IF0 gets the first 858 # value IF1 the second... 859 scan.set_restfreqs(freqs=[1.4e9, 1.67e9]) 1191 # set multiple restfrequencies to all the selected data 1192 scan.set_restfreqs(freqs=[1.4e9, 1.41e9, 1.42e9]) 1193 # If the number of IFs in the data is >= 2 the IF0 gets the first 1194 # value IF1 the second... NOTE that freqs needs to be 1195 # specified in list of list (e.g. [[],[],...] ). 1196 scan.set_restfreqs(freqs=[[1.4e9],[1.67e9]]) 860 1197 #set the given restfrequency for the whole table (by name) 861 1198 scan.set_restfreqs(freqs="OH1667") … … 877 1214 # simple value 878 1215 if isinstance(freqs, int) or isinstance(freqs, float): 879 self._setrestfreqs(freqs, "",unit) 1216 # TT mod 1217 #self._setrestfreqs(freqs, "",unit) 1218 self._setrestfreqs([freqs], [""],unit) 880 1219 # list of values 881 1220 elif isinstance(freqs, list) or isinstance(freqs, tuple): 882 1221 # list values are scalars 883 1222 if isinstance(freqs[-1], int) or isinstance(freqs[-1], float): 1223 self._setrestfreqs(freqs, [""], unit) 1224 # list values are tuples, (value, name) 1225 elif isinstance(freqs[-1], dict): 1226 #sel = selector() 1227 #savesel = self._getselection() 1228 #iflist = self.getifnos() 1229 #for i in xrange(len(freqs)): 1230 # sel.set_ifs(iflist[i]) 1231 # self._setselection(sel) 1232 # self._setrestfreqs(freqs[i], "",unit) 1233 #self._setselection(savesel) 1234 self._setrestfreqs(freqs["value"], 1235 freqs["name"], unit) 1236 elif isinstance(freqs[-1], list) or isinstance(freqs[-1], tuple): 884 1237 sel = selector() 885 1238 savesel = self._getselection() 886 1239 iflist = self.getifnos() 1240 if len(freqs)>len(iflist): 1241 raise ValueError("number of elements in list of list exeeds the current IF selections") 887 1242 for i in xrange(len(freqs)): 888 1243 sel.set_ifs(iflist[i]) 889 1244 self._setselection(sel) 890 self._setrestfreqs(freqs[i], "",unit) 891 self._setselection(savesel) 892 # list values are tuples, (value, name) 893 elif isinstance(freqs[-1], dict): 894 sel = selector() 895 savesel = self._getselection() 896 iflist = self.getifnos() 897 for i in xrange(len(freqs)): 898 sel.set_ifs(iflist[i]) 899 self._setselection(sel) 900 self._setrestfreqs(freqs[i]["value"], 901 freqs[i]["name"], "MHz") 1245 self._setrestfreqs(freqs[i], [""], unit) 902 1246 self._setselection(savesel) 903 1247 # freqs are to be taken from a linecatalog … … 905 1249 sel = selector() 906 1250 savesel = self._getselection() 907 iflist = self.getifnos()908 1251 for i in xrange(freqs.nrow()): 909 1252 sel.set_ifs(iflist[i]) … … 963 1306 msg = "Illegal file name '%s'." % (filename) 964 1307 if rcParams['verbose']: 965 print msg 1308 #print msg 1309 asaplog.push( msg ) 1310 print_log( 'ERROR' ) 966 1311 else: 967 1312 raise IOError(msg) … … 1020 1365 except RuntimeError, msg: 1021 1366 if rcParams['verbose']: 1022 print msg 1367 #print msg 1368 print_log() 1369 asaplog.push( str(msg) ) 1370 print_log( 'ERROR' ) 1023 1371 return 1024 1372 else: raise 1025 1373 s._add_history("average_time", varlist) 1374 print_log() 1026 1375 return s 1027 1376 … … 1051 1400 s = scantable(self._math._convertflux(self, d, eta, jyperk)) 1052 1401 s._add_history("convert_flux", varlist) 1402 print_log() 1053 1403 if insitu: self._assign(s) 1054 1404 else: return s … … 1103 1453 s = scantable(self._math._gainel(self, poly, filename, method)) 1104 1454 s._add_history("gain_el", varlist) 1455 print_log() 1105 1456 if insitu: 1106 1457 self._assign(s) … … 1130 1481 s = scantable(self._math._freq_align(self, reftime, method)) 1131 1482 s._add_history("freq_align", varlist) 1483 print_log() 1132 1484 if insitu: self._assign(s) 1133 1485 else: return s … … 1158 1510 s = scantable(self._math._opacity(self, tau)) 1159 1511 s._add_history("opacity", varlist) 1512 print_log() 1160 1513 if insitu: self._assign(s) 1161 1514 else: return s … … 1176 1529 s = scantable(self._math._bin(self, width)) 1177 1530 s._add_history("bin", varlist) 1531 print_log() 1178 1532 if insitu: 1179 1533 self._assign(s) … … 1200 1554 s = scantable(self._math._resample(self, method, width)) 1201 1555 s._add_history("resample", varlist) 1556 print_log() 1202 1557 if insitu: self._assign(s) 1203 1558 else: return s … … 1218 1573 s = scantable(self._math._averagepol(self, mask, weight.upper())) 1219 1574 s._add_history("average_pol", varlist) 1575 print_log() 1220 1576 return s 1221 1577 … … 1235 1591 s = scantable(self._math._averagebeams(self, mask, weight.upper())) 1236 1592 s._add_history("average_beam", varlist) 1593 print_log() 1237 1594 return s 1238 1595 … … 1263 1620 except RuntimeError, msg: 1264 1621 if rcParams['verbose']: 1265 print msg 1622 #print msg 1623 print_log() 1624 asaplog.push( str(msg) ) 1625 print_log( 'ERROR' ) 1266 1626 return 1267 1627 else: 1268 1628 raise 1269 1629 s._add_history("convert_pol", varlist) 1630 print_log() 1270 1631 return s 1271 1632 1272 1633 @print_log_dec 1273 def smooth(self, kernel="hanning", width=5.0, order=2, insitu=None):1634 def smooth(self, kernel="hanning", width=5.0, order=2, plot=False, insitu=None): 1274 1635 """ 1275 1636 Smooth the spectrum by the specified kernel (conserving flux). … … 1286 1647 specify the order of the polnomial. Ignored by all other 1287 1648 kernels. 1649 plot: plot the original and the smoothed spectra. 1650 In this each indivual fit has to be approved, by 1651 typing 'y' or 'n' 1288 1652 insitu: if False a new scantable is returned. 1289 1653 Otherwise, the scaling is done in-situ … … 1295 1659 self._math._setinsitu(insitu) 1296 1660 varlist = vars() 1661 1662 if plot: orgscan = self.copy() 1663 1297 1664 s = scantable(self._math._smooth(self, kernel.lower(), width, order)) 1298 1665 s._add_history("smooth", varlist) 1666 1667 if plot: 1668 if rcParams['plotter.gui']: 1669 from asap.asaplotgui import asaplotgui as asaplot 1670 else: 1671 from asap.asaplot import asaplot 1672 self._p=asaplot() 1673 self._p.set_panels() 1674 ylab=s._get_ordinate_label() 1675 #self._p.palette(0,["#777777","red"]) 1676 for r in xrange(s.nrow()): 1677 xsm=s._getabcissa(r) 1678 ysm=s._getspectrum(r) 1679 xorg=orgscan._getabcissa(r) 1680 yorg=orgscan._getspectrum(r) 1681 self._p.clear() 1682 self._p.hold() 1683 self._p.set_axes('ylabel',ylab) 1684 self._p.set_axes('xlabel',s._getabcissalabel(r)) 1685 self._p.set_axes('title',s._getsourcename(r)) 1686 self._p.set_line(label='Original',color="#777777") 1687 self._p.plot(xorg,yorg) 1688 self._p.set_line(label='Smoothed',color="red") 1689 self._p.plot(xsm,ysm) 1690 ### Ugly part for legend 1691 for i in [0,1]: 1692 self._p.subplots[0]['lines'].append([self._p.subplots[0]['axes'].lines[i]]) 1693 self._p.release() 1694 ### Ugly part for legend 1695 self._p.subplots[0]['lines']=[] 1696 res = raw_input("Accept smoothing ([y]/n): ") 1697 if res.upper() == 'N': 1698 s._setspectrum(yorg, r) 1699 self._p.unmap() 1700 self._p = None 1701 del orgscan 1702 1703 print_log() 1299 1704 if insitu: self._assign(s) 1300 1705 else: return s … … 1321 1726 """ 1322 1727 if insitu is None: insitu = rcParams['insitu'] 1728 if not insitu: 1729 workscan = self.copy() 1730 else: 1731 workscan = self 1323 1732 varlist = vars() 1324 1733 if mask is None: 1325 1734 mask = [True for i in xrange(self.nchan(-1))] 1735 1326 1736 from asap.asapfitter import fitter 1327 1737 try: 1328 1738 f = fitter() 1329 f.set_scan(self, mask)1330 #f.set_function(poly=order)1331 1739 if uselin: 1332 1740 f.set_function(lpoly=order) 1333 1741 else: 1334 1742 f.set_function(poly=order) 1335 s = f.auto_fit(insitu, plot=plot) 1336 s._add_history("poly_baseline", varlist) 1337 if insitu: self._assign(s) 1338 else: return s 1743 1744 rows = range(workscan.nrow()) 1745 if len(rows) > 0: 1746 self.blpars = [] 1747 1748 for r in rows: 1749 # take into account flagtra info (CAS-1434) 1750 flagtra = workscan._getmask(r) 1751 actualmask = mask[:] 1752 if len(actualmask) == 0: 1753 actualmask = list(flagtra[:]) 1754 else: 1755 if len(actualmask) != len(flagtra): 1756 raise RuntimeError, "Mask and flagtra have different length" 1757 else: 1758 for i in range(0, len(actualmask)): 1759 actualmask[i] = actualmask[i] and flagtra[i] 1760 f.set_scan(workscan, actualmask) 1761 f.x = workscan._getabcissa(r) 1762 f.y = workscan._getspectrum(r) 1763 f.data = None 1764 f.fit() 1765 if plot: 1766 f.plot(residual=True) 1767 x = raw_input("Accept fit ( [y]/n ): ") 1768 if x.upper() == 'N': 1769 self.blpars.append(None) 1770 continue 1771 workscan._setspectrum(f.fitter.getresidual(), r) 1772 self.blpars.append(f.get_parameters()) 1773 1774 if plot: 1775 f._p.unmap() 1776 f._p = None 1777 workscan._add_history("poly_baseline", varlist) 1778 print_log() 1779 if insitu: self._assign(workscan) 1780 else: return workscan 1339 1781 except RuntimeError: 1340 1782 msg = "The fit failed, possibly because it didn't converge." 1341 1783 if rcParams['verbose']: 1342 print msg 1784 #print msg 1785 print_log() 1786 asaplog.push( str(msg) ) 1787 print_log( 'ERROR' ) 1343 1788 return 1344 1789 else: 1345 1790 raise RuntimeError(msg) 1791 1346 1792 1347 1793 def auto_poly_baseline(self, mask=[], edge=(0, 0), order=0, … … 1427 1873 1428 1874 rows = range(workscan.nrow()) 1875 # Save parameters of baseline fits & masklists as a class attribute. 1876 # NOTICE: It does not reflect changes in scantable! 1877 if len(rows) > 0: 1878 self.blpars=[] 1879 self.masklists=[] 1429 1880 asaplog.push("Processing:") 1430 1881 for r in rows: … … 1441 1892 curedge = edge[workscan.getif(r)] 1442 1893 1894 # take into account flagtra info (CAS-1434) 1895 flagtra = workscan._getmask(r) 1896 actualmask = mask[:] 1897 if len(actualmask) == 0: 1898 actualmask = list(flagtra[:]) 1899 else: 1900 if len(actualmask) != len(flagtra): 1901 raise RuntimeError, "Mask and flagtra have different length" 1902 else: 1903 for i in range(0, len(actualmask)): 1904 actualmask[i] = actualmask[i] and flagtra[i] 1905 1443 1906 # setup line finder 1444 fl.find_lines(r, mask, curedge) 1445 f.set_data(workscan._getabcissa(r), workscan._getspectrum(r), 1446 mask_and(workscan._getmask(r), fl.get_mask())) 1907 fl.find_lines(r, actualmask, curedge) 1908 outmask=fl.get_mask() 1909 f.set_scan(workscan, fl.get_mask()) 1910 f.x = workscan._getabcissa(r) 1911 f.y = workscan._getspectrum(r) 1912 f.data = None 1447 1913 f.fit() 1448 x = f.get_parameters() 1914 1915 # Show mask list 1916 masklist=workscan.get_masklist(fl.get_mask(),row=r) 1917 msg = "mask range: "+str(masklist) 1918 asaplog.push(msg, False) 1919 1449 1920 if plot: 1450 1921 f.plot(residual=True) 1451 1922 x = raw_input("Accept fit ( [y]/n ): ") 1452 1923 if x.upper() == 'N': 1924 self.blpars.append(None) 1925 self.masklists.append(None) 1453 1926 continue 1927 1454 1928 workscan._setspectrum(f.fitter.getresidual(), r) 1929 self.blpars.append(f.get_parameters()) 1930 self.masklists.append(masklist) 1455 1931 if plot: 1456 1932 f._p.unmap() … … 1476 1952 self._math._rotate_linpolphase(self, angle) 1477 1953 self._add_history("rotate_linpolphase", varlist) 1954 print_log() 1478 1955 return 1479 1956 … … 1492 1969 self._math._rotate_xyphase(self, angle) 1493 1970 self._add_history("rotate_xyphase", varlist) 1971 print_log() 1494 1972 return 1495 1973 … … 1503 1981 self._math._swap_linears(self) 1504 1982 self._add_history("swap_linears", varlist) 1983 print_log() 1505 1984 return 1506 1985 … … 1513 1992 self._math._invert_phase(self) 1514 1993 self._add_history("invert_phase", varlist) 1994 print_log() 1515 1995 return 1516 1996 … … 1530 2010 s = scantable(self._math._unaryop(self, offset, "ADD", False)) 1531 2011 s._add_history("add", varlist) 2012 print_log() 1532 2013 if insitu: 1533 2014 self._assign(s) … … 1540 2021 Return a scan where all spectra are scaled by the give 'factor' 1541 2022 Parameters: 1542 factor: the scaling factor 2023 factor: the scaling factor (float or 1D float list) 1543 2024 insitu: if False a new scantable is returned. 1544 2025 Otherwise, the scaling is done in-situ … … 1550 2031 self._math._setinsitu(insitu) 1551 2032 varlist = vars() 1552 s = scantable(self._math._unaryop(self, factor, "MUL", tsys)) 2033 s = None 2034 import numpy 2035 if isinstance(factor, list) or isinstance(factor, numpy.ndarray): 2036 if isinstance(factor[0], list) or isinstance(factor[0], numpy.ndarray): 2037 from asapmath import _array2dOp 2038 s = _array2dOp( self.copy(), factor, "MUL", tsys ) 2039 else: 2040 s = scantable( self._math._arrayop( self.copy(), factor, "MUL", tsys ) ) 2041 else: 2042 s = scantable(self._math._unaryop(self.copy(), factor, "MUL", tsys)) 1553 2043 s._add_history("scale", varlist) 2044 print_log() 1554 2045 if insitu: 1555 2046 self._assign(s) … … 1594 2085 1595 2086 @print_log_dec 1596 def auto_quotient(self, preserve=True, mode='paired' ):2087 def auto_quotient(self, preserve=True, mode='paired', verify=False): 1597 2088 """ 1598 2089 This function allows to build quotients automatically. 1599 It assumes the observation to have the same num er of2090 It assumes the observation to have the same number of 1600 2091 "ons" and "offs" 1601 2092 Parameters: … … 1634 2125 s = scantable(self._math._auto_quotient(self, mode, preserve)) 1635 2126 s._add_history("auto_quotient", varlist) 2127 print_log() 1636 2128 return s 1637 2129 … … 1656 2148 q = quotient(on, off, preserve) 1657 2149 q._add_history("mx_quotient", varlist) 2150 print_log() 1658 2151 return q 1659 2152 … … 1674 2167 s = scantable(self._math._freqswitch(self)) 1675 2168 s._add_history("freq_switch", varlist) 2169 print_log() 1676 2170 if insitu: self._assign(s) 1677 2171 else: return s … … 1688 2182 self._recalcazel() 1689 2183 self._add_history("recalc_azel", varlist) 2184 print_log() 1690 2185 return 1691 2186 … … 1765 2260 fit = asapfit(self._getfit(row)) 1766 2261 if rcParams['verbose']: 1767 print fit 2262 #print fit 2263 asaplog.push( '%s' %(fit) ) 2264 print_log() 1768 2265 return 1769 2266 else: … … 1862 2359 return (sum(nchans)/len(nchans) == nchans[0]) 1863 2360 1864 def _fill(self, names, unit, average ):2361 def _fill(self, names, unit, average, getpt, antenna): 1865 2362 import os 1866 2363 from asap._asap import stfiller … … 1874 2371 if rcParams['verbose']: 1875 2372 asaplog.push(msg) 1876 print asaplog.pop().strip() 2373 #print asaplog.pop().strip() 2374 print_log( 'ERROR' ) 1877 2375 return 1878 2376 raise IOError(msg) … … 1889 2387 asaplog.push(msg, False) 1890 2388 print_log() 1891 r._open(name, -1, -1)2389 r._open(name, antenna, -1, -1, getpt) 1892 2390 r._read() 1893 2391 if average: … … 1901 2399 if unit is not None: 1902 2400 self.set_fluxunit(unit) 1903 self.set_freqframe(rcParams['scantable.freqframe'])2401 #self.set_freqframe(rcParams['scantable.freqframe']) 1904 2402 1905 2403 def __getitem__(self, key):
Note: See TracChangeset
for help on using the changeset viewer.