Convenience function for common coordinate functions using pyephem

pyephem_coordinate_transforms.py

"""Coordinate transformation methods using PyEphem

http://rhodesmill.org/pyephem/
"""

import numpy as np
import ephem

def transform(in_position, in_system, out_system, observer=None):
    """Simple coordinate transforms between commonly used coordinate systems.
    
    position: input position as a float tuple (longitude, latitude) in degrees
    in_system, out_system: coordinate system strings; one of
        horizon:    Alt,  Az
        equatorial: RA,   DEC
        galactic:   GLON, GLAT
    observer: ephem.Observer (needed when converting to / from horizon system
    
    Returns: transformed input position as a float tuple (longitude, latitude) in degrees

    See http://stackoverflow.com/questions/11169523/how-to-compute-alt-az-for-given-galactic-coordinate-glon-glat-with-pyephem
    """
    # Set the default observer
    observer = observer if observer else HESS()
    # Internally use radians;
    in_position = np.radians(in_position)
    
    # Transform in_position to Equatorial coordinates ra, dec:
    if in_system == 'horizon':
        ra, dec = map(float, observer.radec_of(in_position[0], in_position[1]))
    elif in_system == 'equatorial':
        ra, dec = in_position
    elif in_system == 'galactic':
        galactic = ephem.Galactic(in_position[0], in_position[1])
        equatorial = ephem.Equatorial(galactic)
        ra, dec = equatorial.ra.real, equatorial.dec.real
    else:
        raise RuntimeError('in_system = %s not supported' % in_system)

    # Here we have ra, dec in radians as floats

    # Now transform Equatorial coordinates to out_system:
    if out_system == 'horizon':
        equatorial = ephem.Equatorial(ra, dec)
        body = ephem.FixedBody()
        body._ra = equatorial.ra
        body._dec = equatorial.dec
        body._epoch = equatorial.epoch
        body.compute(observer)
        out_position = body.az, body.alt
    elif out_system == 'equatorial':
        out_position = ra, dec
    elif out_system == 'galactic':
        equatorial = ephem.Equatorial(ra, dec)
        galactic = ephem.Galactic(equatorial)
        out_position = galactic.lon.real, galactic.lat.real
    else:
        raise RuntimeError('out_system = %s not supported' % out_system)
    
    # Clip longitude to 0 .. 360 deg range
    if out_position[0] > 360:
        out_position[0] = out_position[0] - 360

    # Return out position in degrees
    return np.degrees(out_position)    

def test_transform():
    """Test all coordinate transform methods for one test example.
    
    For three systems there are six possible transformations.
    
    @todo: Rewrite these tests as unit tests
    """
    # Example observer:
    observer = HESS()
    # Example coordinate:
    galactic0 = 0., 0.
    print 'galactic0:', galactic0

    print('\ngalactic <-> equatorial:')
    equatorial1 = transform(galactic0, 'galactic', 'equatorial')
    print 'equatorial1:', equatorial1
    galactic1 = transform(equatorial1, 'equatorial', 'galactic')
    print 'galactic1:', galactic1

    print('\ngalactic <-> horizon:')
    horizon1 = transform(galactic0, 'galactic', 'horizon', observer)
    print 'horizon1:', horizon1
    galactic2 = transform(horizon1, 'horizon', 'galactic', observer)
    print 'galactic2:', galactic2

    print('\nequatorial <-> horizon')
    horizon2 = transform(equatorial1, 'equatorial', 'horizon', observer)
    print 'horizon2:', horizon2
    equatorial2 = transform(horizon2, 'horizon', 'equatorial', observer)
    print 'equatorial2:', equatorial2

if __name__ == '__main__':
    test_transform()