A partial implementation of the astropy coordinates API using pyast

skyposition.py

from .angles import Angle
from .. import units as u
import starlink.Ast as Ast

__all__ = [ 'SkyPosition', 'ICRSPosition', 'AZELPosition',
            'ECLIPTICPosition', 'FK4Position', 'FK4_NO_EPosition',
	    'FK5Position', 'GALACTICPosition', 'GEOCENTRICPosition',
	    'HELIOECLIPTICPosition', 'J2000Position',
	    'SUPERGALACTICPosition', 'UNKNOWNPosition' ]


# -----------------------------------------------------------------
# The base class for all types of celestial coordinate positions.
class SkyPosition(object):
    def __init__( self, *args, **kwargs ):

        # First main task: create an AST SkyFrame with default properties. This
        # records all the metadata associated with the requested coordinate system.
        self._skyframe = Ast.SkyFrame()

        # Override the default SkyFrame properties with the supplied values
        self._set_ast_attribute( "system", kwargs )
        self._set_ast_attribute( "epoch", kwargs )
        self._set_ast_attribute( "equinox", kwargs )
        self._set_ast_attribute( "dut1", kwargs )
        self._set_ast_attribute( "obsalt", kwargs )
        self._set_ast_attribute( "obslat", kwargs )
        self._set_ast_attribute( "obslon", kwargs )
        self._set_ast_attribute( "fmtlon", kwargs, attr="Format(1)" )
        self._set_ast_attribute( "fmtlat", kwargs, attr="Format(2)" )

        # If "system" was not supplied, attempt to determine the system
        # from the names of the supplied arguments (e.g "ra,dec"->ICRS,
        # "slon/slat"->supergalactic, etc). This also copies the corresponding
        # longitude and latitude argument values to the "xlon" and "xlat" entries
        # in kwargs.
        if not self._skyframe.test("System"):
            system = None
            system = self._select_system( kwargs, "AZEL", "az", "el", system )
            system = self._select_system( kwargs, "ECLIPTIC", "elon", "elat", system )
            system = self._select_system( kwargs, "FK4", "ra_fk4", "dec_fk4", system )
            system = self._select_system( kwargs, "FK4_NO_E", "ra_fk4_noe", "dec_fk4_noe", system )
            system = self._select_system( kwargs, "FK5", "ra_fk5", "dec_fk5", system )
            system = self._select_system( kwargs, "GALACTIC", "glon", "glat", system )
            system = self._select_system( kwargs, "GALACTIC", "l", "b", system )
            system = self._select_system( kwargs, "GEOCENTRIC", "ra_geo", "dec_geo", system )
            system = self._select_system( kwargs, "HELIOECLIPTIC", "hlon", "hlat", system )
            system = self._select_system( kwargs, "ICRS", "ra", "dec", system )
            system = self._select_system( kwargs, "J2000", "ra_j2000", "dec_j2000", system )
            system = self._select_system( kwargs, "SUPERGALACTIC", "slon", "slat", system )
            system = self._select_system( kwargs, "UNKNOWN", "lon", "lat", system )

            # If a system was found, store it in the SkyFrame and copy the "xlon"
            # and "xlat" arguments to "lon" and "lat", so that we do not need to
            # differentiate between explicitly and implicitly specified systems below.
            if system != None:
               self._skyframe.System = system
               kwargs[ 'lon' ] = kwargs[ 'xlon' ]
               kwargs[ 'lat' ] = kwargs[ 'xlat' ]

        # Next main task: get the floating point longitude and latitude values,
        # in radians, from the supplied arguments. These are stored in
        # self._lonval and self._latval, and are used when # calling AST
        # functions, which required radians. First, get any supplied units,
        # and ensure they are UnitBase objects rather than strings.
        units = [ None, None ]
        if 'unit' in kwargs:
            unit = kwargs.pop( 'unit' )

            if isinstance(unit, tuple) or isinstance(unit, list):
                units = unit
                if len(unit) > 2:
                    raise ValueError('Cannot give more than 2 units while '
                                     'initializing a coordinate')
            else:
               units[ 0 ] = unit
               units[ 1 ] = unit

            if units[ 0 ] and isinstance( units[ 0 ], str ):
               units[ 0 ] = u.Unit( units[ 0 ] )
            if units[ 1 ] and isinstance( units[ 1 ], str ):
               units[ 1 ] = u.Unit( units[ 0 ] )


        # Indicate we have not yet got a lon/lat pair (in radians)
        self._lonval = None
        self._latval = None

        # First see if a lon/lat pair was supplied as a single string via
        # argument "coordstr", or as the first positional argument. If so, use
        # the SkyFrame to convert these strings to radians.
        coordstr =  kwargs.pop( "coordstr", None )
        if coordstr == None and len(args) == 1 and isinstance( args[ 0 ], str ):
            coordstr = args[ 0 ]

        if coordstr != None:
            if isinstance( coordstr, str ):
                (self._lonval,nc1) = self._radians_from_string( coordstr, 1, units[ 0 ] )
                (self._latval,nc2) = self._radians_from_string( coordstr[nc1:], 2, units[ 1 ] )
                if nc1 == 0 or nc2 == 0:
                    raise ValueError("Could not parse {0}/{1} values from the "
                                     "supplied coordstr string: '{2}'.".
                                      format( self._skyframe.Label_1,
                                              self._skyframe.Label_2,
                                              coordstr ))
            else:
                raise ValueError("Supplied coordstr value has type {0} - must "
                                 "be a string.".format(type(coordstr).__name__))

        # If we have two positional arguments, use them as lon and lat values.
        if len(args) == 2:
            if 'lon' in kwargs:
                raise ValueError("The {0} value was specified by more than "
                                "one argument.".format(self._skyframe.Label_1) )
            else:
                kwargs[ 'lon' ] = args[ 0 ];

            if 'lat' in kwargs:
                raise ValueError("The {0} value was specified by more than "
                                "one argument.".format(self._skyframe.Label_2) )
            else:
                kwargs[ 'lat' ] = args[ 1 ];

        # If a lon value was supplied explicitly via argument "lon", raise an
        # error if we already have a lon value, or convert to radians otherwise.
        lon = kwargs.pop( "lon", None )
        if lon != None:
            if  self._lonval == None:

                self._lonval = self._radians_from_object( lon, 1, units[ 0 ] )
            else:
                raise ValueError("The {0} value was specified by more than "
                                "one argument.".format(self._skyframe.Label_1) )

        # If a lat value was supplied explicitly via argument "lat", raise an
        # error if we already have a lat value, or convert to radians otherwise.
        lat = kwargs.pop( "lat", None )
        if lat != None:
            if  self._latval == None:
               self._latval = self._radians_from_object( lat, 2, units[ 1 ] )
            else:
               raise ValueError("The {0} value was specified by more than "
                                "one argument.".format(self._skyframe.Label_2) )

        #  Check we have both lon and lat now.
        if self._lonval == None:
           raise ValueError("No {0} value was specified.".format(self._skyframe.Label_1) )
        elif self._latval == None:
           raise ValueError("No {0} value was specified.".format(self._skyframe.Label_2) )

        # Next main task: create a dict to hold the position in each supported system,
        # as a pair of Angle objects. Initialise it to hold the position in
        # the supplied system. Other systems will be added as they are requested.
        system = self._skyframe.System
        self._lon_angle = { system : Angle( self._lonval, unit=u.radian ) }
        self._lat_angle = { system : Angle( self._latval, unit=u.radian ) }


    # -------------------------------------------------------------------
    # Return the angular separation between two SkyPositions as an
    # Angle. The supplied SkyPosition is converted into the coordinate
    # system of self before finding the separation.
    def separation( self, other ):

       #  Convert the lon/lat values associated with "other" into the coodinate
       #  system of self. These are in radians.
       ( newlon, newlat ) = other._convert( self._skyframe.System, other._lonval,
                                            other._latval )

       #  Get the angular separation in radians.
       sep = self._skyframe.distance( [self._lonval, self._latval],
                                      [other._lonval, other._latval] )
       # Return an Angle.
       return Angle( sep, unit=u.radian )


    # -------------------------------------------------------------------
    # Returns a new SkyPosition that is offset a specified angular distance
    # (given by "offset" - an Angle), along the geodesic curve away from
    # "self" towards SkyPosition "point".
    def offset( self, point, offset ):

        # Check argument data types:
        if not isinstance( point, SkyPosition ):
            raise ValueError( "Point must be a SkyPosition object" )
        if not isinstance( offset, Angle ):
            raise ValueError( "Offset must be an Angle object" )

        # Get the lon/lat at point in radians, converting them to the coordinate
        # system of self.
        p2 = point._convert( self._skyframe.System, point._lonval, point._latval )

        # Find the lon/lat (in radians) of the returned position.
        (newlon, newlat ) = self._skyframe.offset( (self._lonval,self._latval), p2,
                                                   offset.radians )

        # Return the offset SkyPosition.
        return self._make_new( newlon, newlat )


    # -------------------------------------------------------------------
    # Properties that give access to the SkyPosition longitude and latitude values
    # in any coordinate system.
    @property
    def az( self ):
        return self._make_angles( "AZEL", 1 )

    @property
    def el( self ):
        return self._make_angles( "AZEL", 2 )

    @property
    def elon( self ):
        return self._make_angles( "ECLIPTIC", 1 )

    @property
    def elat( self ):
        return self._make_angles( "ECLIPTIC", 2 )

    @property
    def ra_fk4( self ):
        return self._make_angles( "FK4", 1 )

    @property
    def dec_fk4( self ):
        return self._make_angles( "FK4", 2 )

    @property
    def ra_fk4_noe( self ):
        return self._make_angles( "FK4_NO_E", 1 )

    @property
    def dec_fk4_noe( self ):
        return self._make_angles( "FK4_NO_E", 2 )

    @property
    def ra_fk5( self ):
        return self._make_angles( "FK5", 1 )

    @property
    def dec_fk5( self ):
        return self._make_angles( "FK5", 2 )

    @property
    def glon( self ):
        return self._make_angles( "GALACTIC", 1 )

    @property
    def glat( self ):
        return self._make_angles( "GALACTIC", 2 )

    @property
    def l( self ):
        return self._make_angles( "GALACTIC", 1 )

    @property
    def b( self ):
        return self._make_angles( "GALACTIC", 2 )

    @property
    def ra_geo( self ):
        return self._make_angles( "GEOCENTRIC", 1 )

    @property
    def dec_geo( self ):
        return self._make_angles( "GEOCENTRIC", 2 )

    @property
    def hlon( self ):
        return self._make_angles( "HELIOECLIPTIC", 1 )

    @property
    def hlat( self ):
        return self._make_angles( "HELIOECLIPTIC", 2 )

    @property
    def ra( self ):
        return self._make_angles( "ICRS", 1 )

    @property
    def dec( self ):
        return self._make_angles( "ICRS", 2 )

    @property
    def ra_j2000( self ):
        return self._make_angles( "J2000", 1 )

    @property
    def dec_j2000( self ):
        return self._make_angles( "J2000", 2 )

    @property
    def slon( self ):
        return self._make_angles( "SUPERGALACTIC", 1 )

    @property
    def slat( self ):
        return self._make_angles( "SUPERGALACTIC", 2 )

    @property
    def lon( self ):
        return self._make_angles( "UNKNOWN", 1 )

    @property
    def lat( self ):
        return self._make_angles( "UNKNOWN", 2 )



    # -------------------------------------------------------------------
    # Properties that return a new SkyPosition in any coordinate system.
    @property
    def azel( self ):
        return self._make_skyposition( "AZEL" )

    @property
    def ecliptic( self ):
        return self._make_skyposition( "ECLIPTIC" )

    @property
    def fk4( self ):
        return self._make_skyposition( "FK4" )

    @property
    def fk4_no_e( self ):
        return self._make_skyposition( "FK4_NO_E" )

    @property
    def fk5( self ):
        return self._make_skyposition( "FK5" )

    @property
    def galactic( self ):
        return self._make_skyposition( "GALACTIC" )

    @property
    def geocentric( self ):
        return self._make_skyposition( "GEOCENTRIC" )

    @property
    def helioecliptic( self ):
        return self._make_skyposition( "HELIOECLIPTIC" )

    @property
    def icrs( self ):
        return self._make_skyposition( "ICRS" )

    @property
    def j2000( self ):
        return self._make_skyposition( "J2000" )

    @property
    def supergalactic( self ):
        return self._make_skyposition( "SUPERGALACTIC" )

    @property
    def unknown( self ):
        return self._make_skyposition( "UNKNOWN" )



    # -----------------------------------------------------------------
    # Ensure the SkyPosition contains angles for a specified system, and
    # return either the longitude or latitude angle.
    def _make_angles( self, system, axis ):

        # Only calculate the lon/lat values in the new coordinate system if
        # they have not already been found.
        if not system in self._lon_angle:
            ( newlon, newlat ) = self._convert( system, self._lonval, self._latval )
            self._lon_angle[ system ] = Angle( newlon, unit=u.radian )
            self._lat_angle[ system ] = Angle( newlat, unit=u.radian )

        # Select the lon or lat angle to return.
        if axis == 1:
            result = self._lon_angle[ system ]
        else:
            result = self._lat_angle[ system ]

        return result

    # -----------------------------------------------------------------
    # Return a new SkyPosition that represents the same position as # "self", but in
    # a different coordinate system.
    def _make_skyposition( self, system ):
        ( newlon, newlat ) = self._convert( system, self._lonval, self._latval )
        return self._make_new( newlon, newlat, system=system )

    # -----------------------------------------------------------------
    # Return a new SkyPosition that is a copy of self but with the specified
    # lon/lat values (in radians), and system (which default to the system of self)
    def _make_new( self, lon, lat, system=None ):
        result = SkyPosition.__new__( SkyPosition )
        result._skyframe = self._skyframe.copy()
        if system != None:
           result._skyframe.System = system
        result._lonval = lon
        result._latval = lat
        result._lon_angle = { system : Angle( lon, unit=u.radian ) }
        result._lat_angle = { system : Angle( lat, unit=u.radian ) }
        return result


    # -----------------------------------------------------------------
    # Convert a lon/lat pair (in radians) from the coordinate system of
    # self, to another specified system.
    def _convert( self, system, lon, lat ):

        # Create a copy of the Ast.SkyFrame in which self is defined, and set
        # its System property to the required system.
        sf = self._skyframe.copy()
        sf.System = system

        # Attempt to get a mapping from the coordinate system of self to the
        # requested coordinate system. Raise an exception is this was not possible
        # (e.g. this will be the case if the system is "UNKNOWN").
        map = self._skyframe.convert( sf )
        if map == None:
            raise ValueError("Cannot convert SkyPosition {0} to system '{1}'.".format( self, system ) )

        # Otherwise transform the lon/lat values using the mapping
        else:
           newpos = map.tran( [[self._lonval],[self._latval]] )

        # Return the transformed lon/lat pair, in radians.
        return ( newpos[0][0], newpos[1][0] )


    # -----------------------------------------------------------------
    # See if the argument list supplied to the SkyPosition constructor
    # contains values for arguments with names equal to the names of the
    # longitude and latitude axes of a test coordinate system. If so,
    # return the test system as the function value, and store the axis
    # values as arguments "xlon" and "xlat" in the argument list.
    def _select_system( self, kwargs, testsys, lon_name, lat_name, system ):
        if lon_name in kwargs or lat_name in kwargs:
            if system == None:
                system = testsys

                lon =  kwargs.pop( lon_name, None )
                if lon == None:
                   raise ValueError("No {0} value supplied for {1} position.".format( lon_name, system ) )

                lat =  kwargs.pop( lat_name, None )
                if lat == None:
                   raise ValueError("No {0} value supplied for {1} position.".format( lat_name, system ) )

                kwargs['xlon'] = lon
                kwargs['xlat'] = lat
            else:
                raise ValueError("More than one sky position supplied." )
        return system



    # -----------------------------------------------------------------
    # Read the named value from the supplied dict and use it to set the AST
    # SkyFrame attribute with the same name.
    def _set_ast_attribute( self, name, kwargs, attr=None ):
        value =  kwargs.pop( name, None )
        if value != None:
            if attr == None:
               self._skyframe.set("{0}={1}".format(name,value))
            else:
               self._skyframe.set("{0}={1}".format(attr,value))



    # -----------------------------------------------------------------
    # Convert the supplied object to a radian value.
    def _radians_from_object( self, value, axis, unit ):
        result = None
        if isinstance( value, Angle ):
            result = value.radians
        elif isinstance( value, str ):
            ( result, nchar ) = self._radians_from_string( value, axis, unit )
        else:
            result =  self._to_radians( value, unit )
        return result



    # -----------------------------------------------------------------
    # Read a radian angle from a string, returning the radian value and
    # the number of characters read from the string
    def _radians_from_string( self, value, axis, unit ):
        radians = None
        nchar = 0
        if isinstance( value, str ):

            #  If we have a unit, and the supplied string is a number
            #  convert to a floating point value and use the unit to
            #  convert to radians.
            if unit:
                try:
                    radians = self._to_radians( float( value ), unit )
                    nchar = length( value )
                except:
                    pass

            # If the string was not a single numerical value, it may be a sexagisimal 
            # string, so attempt to parse it using the SkyFrame.
            if radians == None:
                ( nchar, radians ) = self._skyframe.unformat( axis, value )

        return ( radians, nchar )



    # -----------------------------------------------------------------
    # Convert a supplied floating point value to radians.
    def _to_radians( self, value, unit ):
        result = None
        if unit is u.hour:
            result = value*15*Ast.DD2R
        elif unit is u.degree:
            result = value*Ast.DD2R
        elif not unit is u.radian:
            raise ValueError("Bad unit {0} for a celestial axis value.".
                              format(unit) )
        return result



    # -----------------------------------------------------------------
    #  Format a SkyPosition for printing.
    def __str__(self):
        return "<{0} {1}={2} ({3}) {4}={5} ({6})>".format(self._skyframe.System,
                                              self._skyframe.Symbol_1,
                                              self._skyframe.format(1,self._lonval),
                                              self._skyframe.Unit_1,
                                              self._skyframe.Symbol_2,
                                              self._skyframe.format(2,self._latval),
                                              self._skyframe.Unit_2)







# -----------------------------------------------------------------
#  Classes for individual types of celestial coordinate positions.
class ICRSPosition(SkyPosition):
    def __init__( self, *args, **kwargs ):
       kwargs['system'] = "ICRS"
       SkyPosition.__init__( self, *args, **kwargs )

class AZELPosition(SkyPosition):
    def __init__( self, *args, **kwargs ):
       kwargs['system'] = "AZEL"
       SkyPosition.__init__( self, *args, **kwargs )

class ECLIPTICPosition(SkyPosition):
    def __init__( self, *args, **kwargs ):
       kwargs['system'] = "ECLIPTIC"
       SkyPosition.__init__( self, *args, **kwargs )

class FK4Position(SkyPosition):
    def __init__( self, *args, **kwargs ):
       kwargs['system'] = "FK4"
       SkyPosition.__init__( self, *args, **kwargs )

class FK4_NO_EPosition(SkyPosition):
    def __init__( self, *args, **kwargs ):
       kwargs['system'] = "FK4_NO_E"
       SkyPosition.__init__( self, *args, **kwargs )

class FK5Position(SkyPosition):
    def __init__( self, *args, **kwargs ):
       kwargs['system'] = "FK5"
       SkyPosition.__init__( self, *args, **kwargs )

class GALACTICPosition(SkyPosition):
    def __init__( self, *args, **kwargs ):
       kwargs['system'] = "GALACTIC"
       SkyPosition.__init__( self, *args, **kwargs )

class GEOCENTRICPosition(SkyPosition):
    def __init__( self, *args, **kwargs ):
       kwargs['system'] = "GEOCENTRIC"
       SkyPosition.__init__( self, *args, **kwargs )

class HELIOECLIPTICPosition(SkyPosition):
    def __init__( self, *args, **kwargs ):
       kwargs['system'] = "HELIOECLIPTIC"
       SkyPosition.__init__( self, *args, **kwargs )

class J2000Position(SkyPosition):
    def __init__( self, *args, **kwargs ):
       kwargs['system'] = "J2000"
       SkyPosition.__init__( self, *args, **kwargs )

class SUPERGALACTICPosition(SkyPosition):
    def __init__( self, *args, **kwargs ):
       kwargs['system'] = "SUPERGALACTIC"
       SkyPosition.__init__( self, *args, **kwargs )

class UNKNOWNPosition(SkyPosition):
    def __init__( self, *args, **kwargs ):
       kwargs['system'] = "UNKNOWN"
       SkyPosition.__init__( self, *args, **kwargs )

skyposition_test.py

#!/usr/bin/python2

from astropy import units as u
from astropy import coordinates as apc


sc = apc.SkyPosition( ra=10.68458, dec=41.26917, unit=u.degree,
                      fmtlon="d.5", fmtlat="d.5" )
print( sc)

sc = apc.SkyPosition( '00h42m44.3s +41d16m9s', fmtlon="d.5", fmtlat="d.5" )
print( sc)

c = apc.SkyPosition(ra=10.68458, dec=41.26917, unit=u.degree)
print( c.ra.hours )
print( c.dec.radians )
print( c.ra.hms )

print( c.galactic )

print( apc.SkyPosition(ra='12h30m49.42s', dec='+12d23m28.044s') )
print( apc.SkyPosition(l=-76.22237, b=74.49108, unit=u.degree) )
print( apc.SkyPosition(az=45.8, el=42.3, unit=u.degree) )
print( apc.ICRSPosition( 187.70592, 12.39112, unit=u.degree ) )
print( apc.FK4Position(187.07317, 12.66715, unit=u.degree) )
print( apc.GALACTICPosition(-76.22237, 74.49108, unit=u.degree ) )
print( apc.ICRSPosition(54.12412, "-41:08:15.162342", unit=u.degree, fmtlon="d.6", fmtlat="d.6" ))
print( apc.ICRSPosition("3:36:29.7888 -41:08:15.162342", unit=u.hour, fmtlon="d.6", fmtlat="d.6" ))
print( apc.ICRSPosition("54.12412 -41:08:15.162342", fmtlon="d.6", fmtlat="d.6" ))
print( apc.ICRSPosition("14.12412 -41:08:15.162342", fmtlon="d.6", fmtlat="d.6" ))


c1 = apc.ICRSPosition('5h23m34.5s -69d45m22s')
c2 = apc.ICRSPosition('0h52m44.8s -72d49m43s')
sep= c1.separation(c2)  #  Angular distance from c1 to c2
print( sep )

c3 = c1.offset( c2, sep )  #  Find point (c3) which is an angle "sep" away from c1 towards c2. 
print( c3.separation(c2) ) #  c2 and c3 should be coincident so this should display zero !

test = apc.SkyPosition(ra='12h30m49.42s', dec='+12d23m28.044s',
                       obslat="22:19:23.2", obslon="155:19:23.2")
print(test.az)   #  Prints (Az,El) corresponding to 12h30m49.42s +12d23m28.044s

This is a partial implementation of the astropy coordinates API, which is supposed to act just as a "proof-of-concept". To avoid confusion between implementations, I talk about "positions" rather than "coordinates", so the main class is the SkyPosition class. The main thing I have missed out for speed is the Cartesian representation of positions, although this could be added. I've probably missed a few other more minor things as well, but this is just a proof of concept.

My implementation adds some extras to the existing API:

1. Several extra coordinate systems: ecliptic, FK4_NO_E, HelioEcliptic, SuperGalactic, J2000 (that's the dynamical version, not FK5), Geocentric Apparent Equatorial, etc.

2. It adds some properties (obslon,obslat,obsalt) that allow conversion to and from horizon coordinates (az,el) and the other coordinate systems. These are documented in the "Attributes" section of http://starlink.jach.hawaii.edu/docs/sun211.htx/node651.html (see also http://starlink.jach.hawaii.edu/docs/sun211.htx/node677.html for extra metadata associated with celestial coordinate systems in AST). Note, obslon, obslat and obsalt all default to zero.

3. I've added a method called "offset" as an example of the extra functions that the AST Frame and SkyFrame classes provides. This new function returns
   a new SkyPosition that is a given distance away from "self" along a geodesic towards another specified SkyPosition.

This gist also contains a short script that runs some of the examples in the API documentation at http://eteq.github.com/astropy/coordinates/index.html#module-astropy.coordinates

Wolfgang has set it up on his Virtual machine. It has an ipython notebook running on http://moria.astro.utoronto.ca:1111 and you can login with pwd astropyuoft.