<commit_before># Copyright 2014 <NAME> <<EMAIL>> and collaborators.
# Licensed under GNU GPL v2.

""" sdmreader -- functions for reading data and metadata from SDM format files

Functions:

read_bdf -- reads data from binary data format and returns numpy array.
calc_uvw -- parses metadata to calculate uvw coordinates for given scan (or time/direction). returns (u,v,w) tuple. Requires CASA libraries.
read_metadata -- parses metadata of SDM file (xml format) to return tuple with two dictionaries (scaninfo, sourceinfo). Scan info defines BDF location per scan.

BDFData class does the heavy lifting to parse binary data format and return numpy array. Does not yet parse flag table.

Note: baseline order used in the bdf is a bit unusual and different from what is assumed when working with a measurement set.
Order of uvw and axis=1 of data array Pythonically would be [i*nants+j for j in range(nants) for i in range(j)], so [ (1,2), (1,3), (2,3), (1,4), ...].
"""

import numpy as np
import os, mmap, math, string, sdmpy, logging
import cPickle as pickle
import xml.etree.ElementTree as et    # sdmpy can do this part...
from email.feedparser import FeedParser
from email.message import Message

logger = logging.getLogger(__name__)

def read_bdf(sdmpath, scan, nskip=0, readints=0, writebdfpkl=False, bdfdir=None):
    """ Reads given range of integrations from sdm of given scan.
    Uses BDFData object to read.
    readints=0 will read all of bdf (skipping nskip).
    Option to write pkl to store bdf info for faster parse next time.
    """

    assert os.path.exists(sdmpath), 'sdmpath %s does not exist' % sdmpath
    scans, sources = read_metadata(sdmpath, scan, bdfdir=bdfdir)
    assert scans[scan]['bdfstr'], 'bdfstr not defined for scan %d' % scan
    bdffile = scans[scan]['bdfstr']

    assert os.path.exists(bdffile), 'Could not find bdf for scan %d and bdfstr %s.' % (scan, scans[scan]['bdfstr'])

    with open(bdffile, 'r') as fp:
        # define bdfpkldir
        if writebdfpkl:
            bdfpkldir = os.path.join(sdmpath, 'bdfpkls')   # make place for bdfpkls, if needed
            if not os.path.exists(bdfpkldir):
                os.makedirs(bdfpkldir)
        else:
            bdfpkldir = ''

        bdf = BDFData(fp, bdfpkldir=bdfpkldir).parse()
        if readints == 0:
            readints = bdf.n_integrations - nskip

        logger.info('Reading %d ints starting at int %d' % (readints, nskip))
        data = np.empty( (readints, bdf.n_baselines, bdf.n_channels, len(bdf.crosspols)), dtype='complex64', order='C')
        for i in xrange(readints):
            data[i] = bdf.get_data ('crossData.bin', i+nskip)
#            flag[i] = bdf.get_data ('flags.bin', i+nskip)  # need to get auto+cross parsing right to implement this

    return data

def calc_uvw(sdmfile, scan=0, datetime=0, radec=()):
    """ Calculates and returns uvw in meters for a given SDM, time, and pointing direction.
    sdmfile is path to sdm directory that includes "Station.xml" file.
    scan is scan number defined by observatory.
    datetime is time (as string) to calculate uvw (format: '2014/09/03/08:33:04.20')
    radec is (ra,dec) as tuple in units of degrees (format: (180., +45.))
    """

    # set up CASA tools
    try:
        import casautil
    except ImportError:
        try:
            import pwkit.environments.casa.util as casautil
        except ImportError:
            logger.info('Cannot find pwkit/casautil. No calc_uvw possible.')
            return

    me = casautil.tools.measures()
    qa = casautil.tools.quanta()
    logger.debug('Accessing CASA libraries with casautil.')

    assert os.path.exists(os.path.join(sdmfile, 'Station.xml')), 'sdmfile %s has no Station.xml file. Not an SDM?' % sdmfile

    # get scan info
    scans, sources = read_metadata(sdmfile, scan)

    # default is to use scan info
    if (datetime == 0) and (len(radec) == 0):
        assert scan != 0, 'scan must be set when using datetime and radec'   # default scan value not valid

        logger.info('Calculating uvw for first integration of scan %d of source %s' % (scan, scans[scan]['source']))
        datetime = qa.time(qa.quantity(scans[scan]['startmjd'],'d'), form="ymd", prec=8)[0]
        sourcenum = [kk for kk in sources.keys() if sources[kk]['source'] == scans[scan]['source']][0]
        direction = me.direction('J2000', str(np.degrees(sources[sourcenum]['ra']))+'deg', str(np.degrees(sources[sourcenum]['dec']))+'deg')

    # secondary case is when datetime is also given
    elif (datetime != 0) and (len(radec) == 0):
        assert scan != 0, 'scan must be set when using datetime and radec'   # default scan value not valid
        assert '/' in datetime, 'datetime must be in yyyy/mm/dd/hh:mm:ss.sss format'

        logger.info('Calculating uvw at %s for scan %d of source %s' % (datetime, scan, scans[scan]['source']))
        sourcenum = [kk for kk in sources.keys() if sources[kk]['source'] == scans[scan]['source']][0]
        direction = me.direction('J2000', str(np.degrees(sources[sourcenum]['ra']))+'deg', str(np.degrees(sources[sourcenum]['dec']))+'deg')

    else:
        assert '/' in datetime, 'datetime must be in yyyy/mm/dd/hh:mm:ss.sss format'
        assert len(radec) == 2, 'radec must be (ra,dec) tuple in units of degrees'

        logger.info('Calculating uvw at %s in direction %s' % (datetime, direction))
        logger.info('This mode assumes all antennas used.')
        ra = radec[0]; dec = radec[1]
        direction = me.direction('J2000', str(ra)+'deg', str(dec)+'deg')

    # define metadata "frame" for uvw calculation
    sdm = sdmpy.SDM(sdmfile)
    telescopename = sdm['ExecBlock'][0]['telescopeName'].strip()
    logger.debug('Found observatory name %s' % telescopename)

    me.doframe(me.observatory(telescopename))
    me.doframe(me.epoch('utc', datetime))
    me.doframe(direction)

    # read antpos
    if scan != 0:
        configid = [row.configDescriptionId for row in sdm['Main']
                    if scan == int(row.scanNumber)][0]
        antidlist = [row.antennaId for row in sdm['ConfigDescription']
                  if configid == row.configDescriptionId][0].split(' ')
        stationidlist = [ant.stationId
                     for antid in antidlist
                     for ant in sdm['Antenna']
                     if antid == row.antennaId]
    else:
        stationidlist = [ant.stationId for ant in sdm['Antenna']]<commit_msg>tested calc_uvw for scan-dep ant list<commit_after># Copyright 2014 <NAME> <<EMAIL>> and collaborators.
# Licensed under GNU GPL v2.

""" sdmreader -- functions for reading data and metadata from SDM format files

Functions:

read_bdf -- reads data from binary data format and returns numpy array.
calc_uvw -- parses metadata to calculate uvw coordinates for given scan (or time/direction). returns (u,v,w) tuple. Requires CASA libraries.
read_metadata -- parses metadata of SDM file (xml format) to return tuple with two dictionaries (scaninfo, sourceinfo). Scan info defines BDF location per scan.

BDFData class does the heavy lifting to parse binary data format and return numpy array. Does not yet parse flag table.

Note: baseline order used in the bdf is a bit unusual and different from what is assumed when working with a measurement set.
Order of uvw and axis=1 of data array Pythonically would be [i*nants+j for j in range(nants) for i in range(j)], so [ (1,2), (1,3), (2,3), (1,4), ...].
"""

import numpy as np
import os, mmap, math, string, sdmpy, logging
import cPickle as pickle
import xml.etree.ElementTree as et    # sdmpy can do this part...
from email.feedparser import FeedParser
from email.message import Message

logger = logging.getLogger(__name__)

def read_bdf(sdmpath, scan, nskip=0, readints=0, writebdfpkl=False, bdfdir=None):
    """ Reads given range of integrations from sdm of given scan.
    Uses BDFData object to read.
    readints=0 will read all of bdf (skipping nskip).
    Option to write pkl to store bdf info for faster parse next time.
    """

    assert os.path.exists(sdmpath), 'sdmpath %s does not exist' % sdmpath
    scans, sources = read_metadata(sdmpath, scan, bdfdir=bdfdir)
    assert scans[scan]['bdfstr'], 'bdfstr not defined for scan %d' % scan
    bdffile = scans[scan]['bdfstr']

    assert os.path.exists(bdffile), 'Could not find bdf for scan %d and bdfstr %s.' % (scan, scans[scan]['bdfstr'])

    with open(bdffile, 'r') as fp:
        # define bdfpkldir
        if writebdfpkl:
            bdfpkldir = os.path.join(sdmpath, 'bdfpkls')   # make place for bdfpkls, if needed
            if not os.path.exists(bdfpkldir):
                os.makedirs(bdfpkldir)
        else:
            bdfpkldir = ''

        bdf = BDFData(fp, bdfpkldir=bdfpkldir).parse()
        if readints == 0:
            readints = bdf.n_integrations - nskip

        logger.info('Reading %d ints starting at int %d' % (readints, nskip))
        data = np.empty( (readints, bdf.n_baselines, bdf.n_channels, len(bdf.crosspols)), dtype='complex64', order='C')
        for i in xrange(readints):
            data[i] = bdf.get_data ('crossData.bin', i+nskip)
#            flag[i] = bdf.get_data ('flags.bin', i+nskip)  # need to get auto+cross parsing right to implement this

    return data

def calc_uvw(sdmfile, scan=0, datetime=0, radec=()):
    """ Calculates and returns uvw in meters for a given SDM, time, and pointing direction.
    sdmfile is path to sdm directory that includes "Station.xml" file.
    scan is scan number defined by observatory.
    datetime is time (as string) to calculate uvw (format: '2014/09/03/08:33:04.20')
    radec is (ra,dec) as tuple in units of degrees (format: (180., +45.))
    """

    # set up CASA tools
    try:
        import casautil
    except ImportError:
        try:
            import pwkit.environments.casa.util as casautil
        except ImportError:
            logger.info('Cannot find pwkit/casautil. No calc_uvw possible.')
            return

    me = casautil.tools.measures()
    qa = casautil.tools.quanta()
    logger.debug('Accessing CASA libraries with casautil.')

    assert os.path.exists(os.path.join(sdmfile, 'Station.xml')), 'sdmfile %s has no Station.xml file. Not an SDM?' % sdmfile

    # get scan info
    scans, sources = read_metadata(sdmfile, scan)

    # default is to use scan info
    if (datetime == 0) and (len(radec) == 0):
        assert scan != 0, 'scan must be set when using datetime and radec'   # default scan value not valid

        logger.info('Calculating uvw for first integration of scan %d of source %s' % (scan, scans[scan]['source']))
        datetime = qa.time(qa.quantity(scans[scan]['startmjd'],'d'), form="ymd", prec=8)[0]
        sourcenum = [kk for kk in sources.keys() if sources[kk]['source'] == scans[scan]['source']][0]
        direction = me.direction('J2000', str(np.degrees(sources[sourcenum]['ra']))+'deg', str(np.degrees(sources[sourcenum]['dec']))+'deg')

    # secondary case is when datetime is also given
    elif (datetime != 0) and (len(radec) == 0):
        assert scan != 0, 'scan must be set when using datetime and radec'   # default scan value not valid
        assert '/' in datetime, 'datetime must be in yyyy/mm/dd/hh:mm:ss.sss format'

        logger.info('Calculating uvw at %s for scan %d of source %s' % (datetime, scan, scans[scan]['source']))
        sourcenum = [kk for kk in sources.keys() if sources[kk]['source'] == scans[scan]['source']][0]
        direction = me.direction('J2000', str(np.degrees(sources[sourcenum]['ra']))+'deg', str(np.degrees(sources[sourcenum]['dec']))+'deg')

    else:
        assert '/' in datetime, 'datetime must be in yyyy/mm/dd/hh:mm:ss.sss format'
        assert len(radec) == 2, 'radec must be (ra,dec) tuple in units of degrees'

        logger.info('Calculating uvw at %s in direction %s' % (datetime, direction))
        logger.info('This mode assumes all antennas used.')
        ra = radec[0]; dec = radec[1]
        direction = me.direction('J2000', str(ra)+'deg', str(dec)+'deg')

    # define metadata "frame" for uvw calculation
    sdm = sdmpy.SDM(sdmfile)
    telescopename = sdm['ExecBlock'][0]['telescopeName'].strip()
    logger.debug('Found observatory name %s' % telescopename)

    me.doframe(me.observatory(telescopename))
    me.doframe(me.epoch('utc', datetime))
    me.doframe(direction)

    # read antpos
    if scan != 0:
        configid = [row.configDescriptionId for row in sdm['Main'] if scan == int(row.scanNumber)][0]
        antidlist = [row.antennaId for row in sdm['ConfigDescription'] if configid == row.configDescriptionId][0].split(' ')[2:]
        stationidlist = [ant.stationId for antid in antidlist for ant in sdm['Antenna'] if antid == ant.antennaId]
    else:
        stationidlist = [ant.stationId for ant in sdm['Antenna']]
{'fraction_non_alphanumeric': 0.12251268045259461, 'fraction_numerical': 0.01131486539211861, 'mean_word_length': 3.2269129287598943, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 9, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 10}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6575595', 'n_tokens_mistral': 4266, 'n_tokens_neox': 4082, 'n_words': 1355}