Python tables.table函数代码示例

    def _get_ra_and_decl_from_ms(self, measurement_set):
        """
        This function uses pyrap to read the ra and declanation from a 
        measurement set (used by expected_fluxes_in_fov). This is a position 
        in the sky. These values are stored in the field.phase_dir in the first
        row. All exceptions thrown are caught and logged, return None if reading
        failed
        """

        table = None;
        field = None;
        ra_and_decl = None;

        try:
            # open the ms, get the phase direction
            table = pt.table(measurement_set)
            field = pt.table(table.getkeyword("FIELD"))
            ra_and_decl = field.getcell("PHASE_DIR", 0)[0]

        except Exception, exception:
            #catch all exceptions and log
            self.logger.error("Error loading FIELD/PHASE_DIR from "
                              "measurementset {0} : {1}".format(measurement_set,
                                                                str(exception)))
            raise exception

def updateMSmetadata(msfile):
  #Update history to show that this script has modified original data
  tms = pt.table(msfile,readonly=False,ack=False)
  th = pt.table(tms.getkeyword('HISTORY'), readonly=False, ack=False)
  nr=th.nrows()
  th.addrows(1)
  tr=th.row()
  tr.put(nr,{'TIME': quantity('today').get('s').get_value(),
             'OBSERVATION_ID':0,
             'MESSAGE': 'Applied polarization modifications',
             'PRIORITY': 'NORMAL',
             'ORIGIN': '%s: version = %s' % (__file__,__version__),
             'OBJECT_ID':0, 
             'APPLICATION':__file__,
             'CLI_COMMAND':sys.argv,
             'APP_PARAMS': ['']})

  if not options.linear:
     #Change metadata information to be circular feeds
     feed = pt.table(tms.getkeyword('FEED'),readonly=False,ack=False)
     for tpart in feed.iter('ANTENNA_ID'):
      tpart.putcell('POLARIZATION_TYPE',0,['R','L'])

     polariz = pt.table(tms.getkeyword('POLARIZATION'),readonly=False,ack=False)
     polariz.putcell('CORR_TYPE',0,[5,6,7,8])
     tms.close()

def copy_column_to_ms(ms, inputcol, outputcol, ms_from=None):
    """
    Copies one column to another, within an MS file or between two MS files

    Parameters
    ----------
    ms : str
        MS file receiving copy
    inputcol : str
        Column name to copy from
    outputcol : str
        Column name to copy to
    ms_from : str, optional
        MS file to copy from. If None, the column is copied internally

    """
    t = pt.table(ms, readonly=False, ack=False)
    if ms_from is not None:
        tf = pt.table(ms_from, readonly=False, ack=False)
        data = tf.getcol(inputcol)
        desc = tf.getcoldesc(inputcol)
    else:
        data = t.getcol(inputcol)
        desc = t.getcoldesc(inputcol)

    # Add the output column if needed
    if outputcol not in t.colnames():
        desc['name'] = outputcol
        t.addcols(desc)

    t.putcol(outputcol, data)
    t.flush()
    t.close()

def load_and_compare_data_sets(ms1, ms2):
    # open the two datasets
    ms1 = pt.table(ms1)
    ms2 = pt.table(ms2)

    #get the amount of rows in the dataset
    n_row = len(ms1.getcol('DATA'))
    n_complex_vis = 4

    # create a target array with the same length as the datacolumn
    div_array = numpy.zeros((n_row, 1, n_complex_vis), dtype=numpy.complex64)
    ms1_array = ms1.getcol('DATA')
    ms2_array = ms2.getcol('CORRECTED_DATA')

    div_max = 0
    for idx in xrange(n_row):
        for idy  in xrange(n_complex_vis):

            div_value = ms1_array[idx][0][idy] - ms2_array[idx][0][idy]
            if numpy.abs(div_value) > numpy.abs(div_max):
                div_max = div_value

            div_array[idx][0][idy] = div_value
    print "maximum different value between measurement sets: {0}".format(div_max)
    # Use a delta of about float precision
    if div_max > 1e-6:
        print "The measurement sets are contained a different value"
        print "failed delta test!"
        return False

    return True

def splitdataset(dataset, interval, out):
    name = dataset.split("/")[-1]
    print "Splitting {0} by {1} sec intervals...".format(name, interval)
    t = pt.table(dataset, ack=False)
    starttime = t[0]["TIME"]
    endtime = t[t.nrows() - 1]["TIME"]
    numberofsplits = int((endtime - starttime) / interval)
    for split in range(0, numberofsplits):
        outputname = os.path.join(out, "splitMS", name + ".{0}sec_{1:04d}.split".format(int(interval), split + 1))
        if split == 0:
            thisstart = starttime - 2.0
        else:
            thisstart = starttime + (float(split) * interval)
        thisend = starttime + ((float(split) + 1) * interval)
        t1 = t.query(
            "TIME > "
            + str(thisstart)
            + " && \
		TIME < "
            + str(thisend),
            sortlist="TIME,ANTENNA1,ANTENNA2",
        )
        t1.copy(outputname, True)
        t1.close()
        if split == 0:
            thisstart += 2.0
        t1 = pt.table(outputname + "/OBSERVATION", ack=False, readonly=False)
        thistimerange = np.array([thisstart, thisend])
        t1.putcell("TIME_RANGE", 0, thistimerange)
        t1.putcell("LOFAR_OBSERVATION_START", 0, thisstart)
        t1.putcell("LOFAR_OBSERVATION_END", 0, thisend)
        t1.close()
    t.close()

def main(inms='',beVerbose=False):
        didWarn = False
        whatSkipped = {}
        t = pt.table(inms, ack=False)
        th = pt.table(t.getkeyword('HISTORY'), ack=False)
        colnames = th.colnames()
        nrows = th.nrows()
        print 'The HISTORY table in %s has %d rows' % (inms, nrows)
        for row in th:
                if row['APPLICATION'] == 'imager' or row['APPLICATION'] == 'OLAP' or row['APPLICATION'] == 'ms':
                        if beVerbose:
                                print '%s was run at time %f with parameters:' % (row['APPLICATION'],row['TIME'])
                                for r in row['APP_PARAMS']:
                                        print '\t%s' % (r)
                        else:
                                if not didWarn:
                                        print '(Skipping OLAP, imager, and ms rows, use -v to print them)'
                                        didWarn = True
                                if row['APPLICATION'] in whatSkipped:
                                        whatSkipped[row['APPLICATION']] += 1
                                else:
                                        whatSkipped[row['APPLICATION']] = 1
                else:
                        print '%s was run at time %f with parameters:' % (row['APPLICATION'],row['TIME'])
                        for r in row['APP_PARAMS']:
                                print '\t%s' % (r)
        print 'Overview of skipped rows:'
        for key in whatSkipped:
                print '\t%s:\tskipped %d times' % (key,whatSkipped[key])

def copy_column_to_bands(mslist, ms_from, inputcol, outputcol):
    """
    Copies one column from an MS file to multiple MS files (bands)

    Parameters
    ----------
    mslist : list
        MS files receiving copy
    ms_from : str
        MS file to copy from.
    inputcol : str
        Column name to copy from
    outputcol : str
        Column name to copy to

    """
    datain = pt.table(ms_from)
    data = datain.getcol(inputcol, nrow=1)
    numberofchans = numpy.int(numpy.shape(data)[1])
    chanperms = numberofchans/numpy.int(len(mslist))

    for ms_id, ms in enumerate(mslist):
        if os.path.isdir(ms):
            data = datain.getcolslice(inputcol, [chanperms*ms_id,0], [(chanperms*(ms_id+1))-1,3])
            dataout = pt.table(ms, readonly=False)
            dataout.putcol(outputcol, data)
            dataout.flush()
            dataout.close()

def gain2matlab(msname='test1.MS', gainfilename='bbsgain.mat', timeslot=0, instrumentname='instrument'):
  parmdb=msname+'/'+instrumentname
  antenna=msname+'/ANTENNA'

  valstable=table(parmdb,ack=False)
  namestable=table(parmdb+"::NAMES",ack=False)
  antennatable=table(antenna,ack=False)

  vals=valstable.col('VALUES')
  names=namestable.col('NAME')
  antennas=antennatable.col('NAME')
  antennamap={};
  for i in range(antennas.nrows()):
    antennamap[antennas[i]]=i

  g = np.zeros((len(antennamap)*2,2),dtype=np.complex)

  for i in range(vals.nrows()):
    (bla, xcor, ycor, reim, ant) = names[i].split(':')
    antnr=antennamap[ant]
    (xcor,ycor)=(int(xcor),int(ycor))
    if reim=="Real":
      val=vals[i][timeslot][0]
    elif reim=="Imag":
      val=vals[i][timeslot][0]*(1.j)
    elif reim=="Phase":
      val=cmath.rect(1,vals[i][timeslot][0])
    #print antnr, xcor, ycor, antnr*2+xcor, ycor, val
    g[antnr*2+ycor][xcor]+=val.conjugate()

  scipy.io.savemat(gainfilename, dict(g=g), oned_as='row')
  print "Stored timeslot",timeslot,"gains from", msname, "/",instrumentname,"as",gainfilename

 def __init__(self,mslist,mss=None):
     """
     mslist is the MS list filename
     """
     import pyrap.tables as pt
     if mss is not None:
         self.mss=mss
         self.mslist=None
     else:
         self.mslist=mslist
         self.mss=[s.strip() for s in open(mslist).readlines()]
     self.obsids = [os.path.basename(ms).split('_')[0] for ms in self.mss]
     self.freqs=[]
     self.channels=[]
     self.hascorrected=[]
     self.dysco=[]
     for ms in self.mss:
         t = pt.table(ms,readonly=True,ack=False)
         colname='CORRECTED_DATA'
         try:
             dummy=t.getcoldesc(colname)
         except RuntimeError:
             dummy=None
         self.hascorrected.append(not(dummy is None))
         self.dysco.append('Dysco' in t.showstructure())
         t.close()
         t = pt.table(ms+'/SPECTRAL_WINDOW', readonly=True, ack=False)
         self.freqs.append(t[0]['REF_FREQUENCY'])
         self.channels.append(t[0]['CHAN_FREQ'])

def read_ms(infile, verbosity=1):
    """ Convert MS to a HDF file
    :param infile:  Measurement Set path
    :return: HDU version of Measurement Set
    """
    pp = PrintLog(verbosity=verbosity)
    ms = pt.table(infile)

    # Create a HDU List for storing HDUs
    hdul = IdiHdulist(verbosity=verbosity)

    # Add each column to the main HDU
    hdu_main = table2hdu(ms, "MAIN", verbosity=verbosity, close_after=False)
    hdul["MAIN"] = hdu_main

    # Now look for other keyword tables
    for key, val in ms.getkeywords().items():
        pp.debug(val)
        if type(val) in (unicode, str):
            if val.startswith("Table: "):
                tblpath = val.strip().split("Table: ")[1]
                pp.h2("Opening %s" % key)
                t = pt.table(tblpath)
                t_hdu = table2hdu(t, key, verbosity=verbosity)
                hdul[key] = t_hdu
        else:
            hdul["MAIN"].header.vals[key] = val

    ms.close()
    return hdul

 def get_summary(self):
     subtables = self.ms.keywordnames()
     for subtable in ('POLARIZATION', 'OBSERVATION', 'FIELD', 
                      'SPECTRAL_WINDOW'):
         if subtable not in subtables:
             sys.stderr.write("Subtable %s missing from MS\n" % subtable)
             sys.exit()
     frequencies = self.get_frequencies()
     self.get_antenntas()
     polarization = {'count': pt.table(
             os.path.join(self.filename, 'POLARIZATION'), ack=False).getcol(
             'NUM_CORR')[0]}
     times = {'time': pt.table(
             os.path.join(self.filename, 'OBSERVATION'), ack=False).getcol(
             'TIME_RANGE')[0]}
     fieldnames = {'fieldnames': pt.table(
             os.path.join(self.filename, 'FIELD'), ack=False).getcol('NAME')}
     phases = {'direction': pt.table(
             os.path.join(self.filename, 'FIELD'), ack=False).getcol('PHASE_DIR')}
     self.summary = {
         'frequencies': frequencies,
         'polarization': polarization,
         'fieldnames': fieldnames,
         'times': times,
         'phases': phases
         }

 def __init__(self, MSfile, timecorr, block, solint, ionfactor, ncores,
     resume, parset, skymodel, parmdb, clobber, solver):
     self.file = MSfile
     self.msname = self.file.split('/')[-1]
     sw = pt.table(self.file + '/SPECTRAL_WINDOW', ack=False)
     self.freq = sw.col('REF_FREQUENCY')[0]
     sw.close()
     obs = pt.table(self.file + '/FIELD', ack=False)
     self.ra = np.degrees(float(obs.col('REFERENCE_DIR')[0][0][0]))
     if self.ra < 0.:
         self.ra=360.+(self.ra)
     self.dec = np.degrees(float(obs.col('REFERENCE_DIR')[0][0][1]))
     obs.close()
     ant = pt.table(self.file + '/ANTENNA', ack=False)
     diam = float(ant.col('DISH_DIAMETER')[0])
     ant.close()
     self.fwhm_deg = 1.1*((3.0e8/self.freq)/diam)*180./np.pi
     self.name = str(self.freq)
     self.timecorr = timecorr
     self.sol_block = block
     self.ionfactor = ionfactor
     self.ncores = ncores
     self.resume = resume
     self.solint = solint
     self.parset = parset
     self.input_parmdb = parmdb
     self.output_parmdb = 'instrument'
     self.skymodel = skymodel
     self.clobber = clobber
     self.solver = solver

    def __init__(self, ms):
        self.timepara={'start':0, 'end':0, 'step':0, 'cent':0}
	self.freqpara={'start':0, 'end':0, 'step':0, 'cent':0}
	self.msname = ms
	if not os.path.isdir(ms): sys.exit('INPUT MS DOES NOT EXIST!')
        ##########Getting Time parameters first#############
	t = pt.table(ms, readonly=True, ack=False)
	t1 = t.sort ('unique desc TIME')
	self.timepara['step'] = t.getcell('EXPOSURE',0)
	self.timepara['start'] =  np.min(t.getcol('TIME'))-self.timepara['step']/2.
	self.timepara['end'] =  np.max(t.getcol('TIME'))+self.timepara['step']/2.
	self.timepara['cent'] = self.timepara['start']+(self.timepara['end']-self.timepara['start'])/2.
	self.mstimevalues = t1.getcol('TIME')[::-1]
	t1.close()
        ##########Getting Frequency Parameters###################
	freq=pt.table(t.getkeyword("SPECTRAL_WINDOW"), readonly=True, ack=False)
	self.fullband = freq.getcell('TOTAL_BANDWIDTH', 0)
	self.freqpara['cent'] = freq.getcell('REF_FREQUENCY', 0)
	self.freqpara['step'] = freq.getcell('CHAN_WIDTH', 0)[0]
        self.msfreqvalues = freq.getcell('CHAN_FREQ', 0)
	self.freqpara['start'] = self.msfreqvalues[0]-self.freqpara['step']/2.
	self.freqpara['end'] = self.msfreqvalues[-1]+self.freqpara['step']/2.
	freq.close()
        ##########Getting Station Names###################
        antennas = pt.table(t.getkeyword("ANTENNA"), readonly=True, ack=False)
        self.stations = antennas.getcol('NAME')
        antennas.close()
	t.close()

def read_ms(logger, msname, ateam, diameter=None):
    def get_station_diameter(table):
        histable = pt.table(table.getkeyword('HISTORY'), ack=False)
        for line in histable.getcell('APP_PARAMS', 0):
            try:
                key, value = line.split("=")
            except:
                pass
            if key == "Observation.antennaSet":
                antenna_set = value
                break
        if antenna_set == "LBA_INNER":
            logger.debug("LBA_INNER mode")
            return STATION_DIAMETER["LBA_INNER"]
        elif antenna_set[:3] == "LBA":
            logger.debug("LBA_(OUTER,SPARSE,X,Y) mode")
            return STATION_DIAMETER["LBA"]
        elif antenna_set[:3] == "HBA":
            logger.debug("HBA mode")
            return STATION_DIAMETER["HBA"]
        else:
            logger.error("Failed to identify antenna set")

    def field_size_ateam(table):
        logging.debug('Computing field size for A-team')
        fieldtable = table.getkeyword('FIELD').split()[1]
        taqloutput = pt.taql("calc from %s calc max(angdist (DELAY_DIR[0,], [%s]))" % (fieldtable, ", ".join(",".join(src) for src in ATEAM))  )
        return taqloutput[0]

    def field_size_nominal(table, wavelength, diameter):
        if not diameter:
            diameter = get_station_diameter(table)
        logger.debug("Station diameter %f m" % diameter)
        return 1.22*wavelength/diameter

    t = pt.table(msname, readonly=True, ack=False)
    interval = t.getcell('INTERVAL', 0)
    swtable = t.getkeyword('SPECTRAL_WINDOW')
    tsw = pt.table(swtable, readonly=True, ack=False)
    freq = tsw.getcell('REF_FREQUENCY', 0)
    wavelength = 299792458./freq
    maxbl = pt.taql("calc sqrt(max([select sumsqr(UVW[0:1]) from %s]))" % msname)[0] / wavelength
    chwidth = tsw.getcell('CHAN_WIDTH', 0)[0]

    if ateam:
        fieldsize = field_size_ateam(t)
    else:
        fieldsize = field_size_nominal(t, wavelength, diameter)

    logger.debug('Frequency is %f MHz'%(freq/1.e6))
    logger.debug('Wavelength is %f m'%(wavelength))
    logger.debug('Maximum baseline length is %f m = %f lambdas'%(maxbl*wavelength,maxbl))
    logger.debug('Integration time is %f sec'%(interval))
    logger.debug('Channel width is %f Hz'%(chwidth))
    logger.debug('Field size is %f degrees'%(fieldsize*180./3.14159))

    return fieldsize, maxbl, freq, interval, chwidth

def rename2(SB, obsid):
	SBtable=pt.table("{0}/OBSERVATION".format(SB), ack=False)
	beam=int(SBtable.col("LOFAR_SUB_ARRAY_POINTING")[0])
	SBtable.close()
	SBtable=pt.table("{0}/SPECTRAL_WINDOW".format(SB), ack=False)
	sbno=int(SBtable.col("NAME")[0].split("-")[-1])
	SBtable.close()
	newname="{0}_SAP{1:03d}_SB{2:03d}_uv.MS.dppp".format(obsid,beam,sbno)
	return newname