ledatelescope / bifrost

A stream processing framework for high-throughput applications.
BSD 3-Clause "New" or "Revised" License
66 stars 29 forks source link

sigproc writer block needs some love #110

Open telegraphic opened 6 years ago

telegraphic commented 6 years ago

The current sigproc reader block a) doesn't treat angles correctly, and b) will crash in most cases! My fix so far has been to use astropy, but that adds a dependency, so starting an issue to discuss merits.

Angles

The sigproc definition for angles is actually a float, HHMMSS.ss. For example, 17h45m00.12s needs to be written as a float 174500.12. This code is ugly but parses hms or dms:

from astropy.coordinates import Angle

def to_sigproc_angle(angle_val):
    """ Convert an astropy.Angle to the ridiculous sigproc angle format string. """
    x = str(angle_val)

    if '.' in x:
        if 'h' in x:
                d, m, s, ss = int(x[0:x.index('h')]), int(x[x.index('h')+1:x.index('m')]), \
                int(x[x.index('m')+1:x.index('.')]), float(x[x.index('.'):x.index('s')])
        if 'd' in x:
            d, m, s, ss = int(x[0:x.index('d')]), int(x[x.index('d')+1:x.index('m')]), \
            int(x[x.index('m')+1:x.index('.')]), float(x[x.index('.'):x.index('s')])
    else:
        if 'h' in x:
            d, m, s = int(x[0:x.index('h')]), int(x[x.index('h')+1:x.index('m')]), \
            int(x[x.index('m')+1:x.index('s')])
        if 'd' in x:
            d, m, s = int(x[0:x.index('d')]), int(x[x.index('d')+1:x.index('m')]), \
            int(x[x.index('m')+1:x.index('s')])
        ss = 0
    num = str(d).zfill(2) + str(m).zfill(2) + str(s).zfill(2)+ '.' + str(ss).split(".")[-1]
    return np.float64(num)

...

# Get RA and DEC and convert to sigproc format
ra = Angle(ihdr['raj'], unit='hourangle')
dec = Angle(ihdr['dej'], unit='degree')
sigproc_hdr['src_raj'] =  to_sigproc_angle(ra)
sigproc_hdr['src_dej'] =  to_sigproc_angle(dec)

Crash bug

The code fails on L204 and L291 as axnames is of type list, while the RHS a tuple.

Easiest fix is:

 if ndim >= 3 and axnames[-3:] == ('time', 'pol', 'freq'):
becomes
 if ndim >= 3 and tuple(axnames[-3:]) == ('time', 'pol', 'freq'):
benbarsdell commented 6 years ago

Thanks for finding these bugs.

Couldn't you use http://docs.astropy.org/en/stable/api/astropy.coordinates.Angle.html#astropy.coordinates.Angle.dms

Astropy does seem like a heavyweight solution for this. How/where do you need to interact with the angle value in the headers? E.g., are you setting them in a source block, loading them from a Sigproc read block, something else? I guess we need to decide on a standard representation of them in headers.

telegraphic commented 6 years ago

Yes, generally read from guppi raw files, which are often strings (unsure what the spec is, or if people actually follow it).

An yes, that code is particularly horrible (not written by me in my defense!), and came from a non-astropy-using codebase.

Astropy angles and units (http://docs.astropy.org/en/v0.2.1/units/index.html) are both pretty useful. I think we either use them extensively (instead of pint?), or not at all...

jaycedowell commented 2 years ago

@telegraphic is this still a thing? It looks like the crash has been fixed but nothing was done about the angle representation.