Coordinates API, take 2

coordinates_api_2.py

@@ -0,0 +1,367 @@
+# -*- coding: utf-8 -*-
+
+from astropy import coordinates as coords
+from astropy import units as u
+
+#<-----------------Classes for representation of coordinate data--------------->
+# These classes inherit from a common base class and internally contain Quantity
+# objects, which are arrays (although they may act as scalars, like numpy's
+# length-0  "arrays")
+
+# They can be initialized with a variety of ways that make intuitive sense.
+# Distance is optional.
+coords.SphericalRepresentation(lon=8*u.hour, lat=5*u.deg)
+coords.SphericalRepresentation(lon=8*u.hourangle, lat=5*u.deg)
+coords.SphericalRepresentation(lon=8*u.hourangle, lat=5*u.deg, distance=10*u.kpc)
+
+# In the initial implementation, the lat/lon/distance arguments to the
+# initializer must be in order. A *possible* future change will be to allow
+# smarter guessing of the order.  E.g. `Latitude` and `Longitude` objects can be
+# given in any order.
+coords.SphericalRepresentation(coords.Longitude(8, u.hour), coords.Latitude(5, u.deg))
+coords.SphericalRepresentation(coords.Longitude(8, u.hour), coords.Latitude(5, u.deg), coords.Distance(10, u.kpc))
+
+# Arrays of any of the inputs are fine
+coords.SphericalRepresentation(lon=[8, 9]*u.hourangle, lat=[5, 6]*u.deg)
+
+# Default is to copy arrays, but optionally, it can be a reference
+coords.SphericalRepresentation(lon=[8, 9]*u.hourangle, lat=[5, 6]*u.deg, copy=False)
+
+# strings are parsed by `Latitude` and `Longitude` constructors, so no need to
+# implement parsing in the Representation classes
+coords.SphericalRepresentation(lon='2h6m3.3s', lat='5rad')
+
+# Or, you can give `Quantity`s with keywords, and they will be internally
+# converted to Angle/Distance
+c1 = coords.SphericalRepresentation(lon=8*u.hourangle, lat=5*u.deg, distance=10*u.kpc)
+
+# Can also give another representation object with the `reprobj` keyword.
+c2 = coords.SphericalRepresentation(reprobj=c1)
+# lat/lon/distance must be None in the case below because it's ambiguous if they
+# should come from the `c1` object or the explicitly-passed keywords.
+c2 = coords.SphericalRepresentation(lon=8*u.hourangle, lat=5*u.deg, reprobj=c1) #raises ValueError
+
+#  distance, lat, and lon typically will just match in shape
+coords.SphericalRepresentation(lon=[8, 9]*u.hourangle, lat=[5, 6]*u.deg, distance=[10, 11]*u.kpc)
+# if the inputs are not the same, if possible they will be broadcast following
+# numpy's standard broadcasting rules.
+c2 = coords.SphericalRepresentation(lon=[8, 9]*u.hourangle, lat=[5, 6]*u.deg, distance=10*u.kpc)
+assert len(c2.distance) == 2
+#when they can't be broadcast, it is a ValueError (same as Numpy)
+with raises(ValueError):
+    c2 = coords.SphericalRepresentation(lon=[8, 9, 10]*u.hourangle, lat=[5, 6]*u.deg)
+
+# It's also possible to pass in scalar quantity lists with mixed units. These
+# are converted to array quantities following the same rule as `Quantity`: all
+# elements are converted to match the first element's units.
+c2 = coords.SphericalRepresentation(lon=[8*u.hourangle, 135*u.deg],
+                                    lat=[5*u.deg, (6*pi/180)*u.rad])
+assert c2.lat.unit == u.deg and c2.lon.unit == u.hourangle
+assert c2.lon[1].value == 9
+
+# The Quantity initializer itself can also be used to force the unit even if the
+# first element doesn't have the right unit
+lon = u.Quantity([120*u.deg, 135*u.deg], u.hourangle)
+lat = u.Quantity([(5*pi/180)*u.rad, 0.4*u.hourangle], u.deg)
+c2 = coords.SphericalRepresentation(lon, lat)
+
+
+# regardless of how input, the `lat` and `lon` come out as angle/distance
+assert isinstance(c1.lat, coords.Angle)
+assert isinstance(c1.lat, coords.Latitude)  # `Latitude` is an `Angle` subclass
+assert isinstance(c1.distance, coords.Distance)
+
+# but they are read-only, as representations are immutable once created
+with raises(AttributeError):
+    c1.lat = coords.Latitude(...)
+# Note that it is still possible to modify the array in-place, but this is not
+# sanctioned by the API, as this would prevent things like caching.
+c1.lat[:] = [0] * u.deg  # possible, but NOT SUPPORTED
+
+# To address the fact that there are various other conventions for how spherical
+# coordinates are defined, other conventions can be included as new classes.
+# Later there may be other conventions that we implement - for now just the
+# physics convention, as it is one of the most common cases.
+c3 = coords.PhysicistSphericalRepresentation(phi=120*u.deg, theta=85*u.deg, rho=3*u.kpc)
+c3 = coords.PhysicistSphericalRepresentation(phi=120*u.deg, theta=85*u.deg, r=3*u.kpc)
+with raises(ValueError):
+    c3 = coords.PhysicistSphericalRepresentation(phi=120*u.deg, theta=85*u.deg, r=3*u.kpc, rho=4*u.kpc)
+
+# first dimension must be length-3 if a lone `Quantity` is passed in.
+c1 = coords.CartesianRepresentation(randn(3, 100) * u.kpc)
+assert c1.xyz.shape[0] == 0
+assert c1.unit == u.kpc
+# using `c1.xyz.unit` is the same as `c1.unit`, so it's not necessary to do this,
+# but it illustrates that `xyz` is a full-fledged `Quantity`.
+assert c1.xyz.unit == u.kpc
+assert c1.x.shape[0] == 100
+assert c1.y.shape[0] == 100
+assert c1.z.shape[0] == 100
+# can also give each as separate keywords
+coords.CartesianRepresentation(x=randn(100)*u.kpc, y=randn(100)*u.kpc, z=randn(100)*u.kpc)
+# if the units don't match but are all distances, they will automatically be
+# converted to match `x`
+xarr, yarr, zarr = randn(3, 100)
+c1 = coords.CartesianRepresentation(x=xarr*u.kpc, y=yarr*u.kpc, z=zarr*u.kpc)
+c2 = coords.CartesianRepresentation(x=xarr*u.kpc, y=yarr*u.kpc, z=zarr*u.pc)
+assert c2.unit == u.kpc
+assert c1.z.kpc / 1000 == c2.z.kpc
+# although if they are not `Distance` compatible, it's an error
+c2 = coords.CartesianRepresentation(x=randn(100)*u.kpc, y=randn(100)*u.kpc, z=randn(100)*u.deg) # raises UnitsError
+
+# CartesianRepresentation can also accept raw arrays and a `unit` keyword
+# instead of having units attached to each of `x`, `y`, and `z`. Note that this
+# is *not* the case for other representations - it's only sensible for
+# Cartesian, because all of the data axes all have the same unit.
+coords.CartesianRepresentation(x=randn(100), y=randn(100), z=randn(100), unit=u.kpc)
+
+# representations convert into other representations via  `represent_as`
+srep = coords.SphericalRepresentation(lon=90*u.deg, lat=0*u.deg, distance=1*u.pc)
+crep = srep.represent_as(coords.CartesianRepresentation)
+assert crep.x.value == 0 and crep.y.value == 1 and crep.z.value == 0
+# The functions that actually do the conversion are defined via methods on the
+# representation classes. This may later be expanded into a full registerable
+# transform graph like the coordinate frames, but initially it will be a simpler
+# method system
+
+
+# Future extensions: add additional useful representations like cylindrical or elliptical
+
+
+#<---------------------Reference Frame/"Low-level" classes--------------------->
+# The low-level classes have a dual role: they act as specifiers of coordinate
+# frames and they *may* also contain data as one of the representation objects,
+# in which case they are the actual coordinate objects themselves.
+
+
+# They can always accept a representation as a first argument
+icrs = coords.ICRS(coords.SphericalRepresentation(lon=8*u.hour, lat=5*u.deg))
+
+# which is stored as the `data` attribute
+assert icrs.data.lat == 5*u.deg
+assert icrs.data.lon == 8*u.hour
+
+# Frames that require additional information like equinoxs or obstimes get them
+# as keyword parameters to the frame constructor.  Where sensible, defaults are
+# used. E.g., FK5 is almost always J2000 equinox
+fk5 = coords.FK5(coords.SphericalRepresentation(lon=8*u.hour, lat=5*u.deg))
+J2000 = astropy.time.Time('J2000',scale='utc')
+fk5_2000 = coords.FK5(coords.SphericalRepresentation(lon=8*u.hour, lat=5*u.deg), equinox=J2000)
+assert fk5.equinox == fk5_2000.equionx
+
+# the information required to specify the frame is immutable
+fk5.equinox = J2001  # raises AttributeError
+
+# As is the representation data.
+with raises(AttributeError):
+    fk5.data = ...
+
+# There is also a class-level attribute that lists the attributes needed to
+# identify the frame.  These include attributes like `equinox` shown above.
+assert FK5.frame_attr_names == ('equinox', 'obstime')  # defined on the *class*
+assert fk5.frame_attr_names == ('equinox', 'obstime')  # and hence also in the objects
+# `frame_attr_names` will mainly be used by the high-level class (discussed
+# below) to allow round-tripping between various frames.  It is also part of the
+# public API for other similar developer / advanced users' use.
+
+# The actual position information is accessed via the representation objects
+assert icrs.represent_as(coords.SphericalRepresentation).lat == 5*u.deg
+assert icrs.spherical.lat == 5*u.deg  #  shorthand for the above
+assert icrs.cartesian.z.value > 0
+
+# Many frames have a "preferred" representation, the one in which they are
+# conventionally described, often with a special name for some of the
+# coordinates. E.g., most equatorial coordinate systems are spherical with RA and
+# Dec. This works simply as a shorthand for the longer form above
+
+assert icrs.ra == 5*u.deg
+assert fk5.dec == 8*u.hour
+
+assert icrs.preferred_representation == coords.SphericalRepresentation
+
+# low-level classes can also be initialized with the preferred names:
+icrs_2 = coords.ICRS(ra=8*u.hour, dec=5*u.deg, distance=1*u.kpc)
+assert icrs == icrs2
+
+# and these are taken as the default if keywords are not given:
+icrs_nokwarg = coords.ICRS(8*u.hour, 5*u.deg, distance=1*u.kpc)
+assert icrs_nokwarg.ra == icrs_2.ra and icrs_nokwarg.dec == icrs_2.dec
+
+# they also are capable of computing on-sky or 3d separations from each other,
+# which will be a direct port of the existing methods:
+coo1 = coords.ICRS(ra=0*u.hour, dec=0*u.deg)
+coo2 = coords.ICRS(ra=0*u.hour, dec=1*u.deg)
+# `separation` is the on-sky separation
+assert coo1.separation(coo2).degree == 1.0
+
+# while `separation_3d` includes the 3D distance information
+coo3 = coords.ICRS(ra=0*u.hour, dec=0*u.deg, distance=1*u.kpc)
+coo4  = coords.ICRS(ra=0*u.hour, dec=0*u.deg, distance=2*u.kpc)
+assert coo3.separation_3d(coo4).kpc == 1.0
+
+# The next example fails because `coo1` and `coo2` don't have distances
+assert coo1.separation_3d(coo2).kpc == 1.0  # raises ValueError
+
+
+# the frames also know how to give a reasonable-looking string of themselves,
+# based on the preferred representation and possibly distance
+assert str(icrs_2) == '<ICRS RA=120.000 deg, Dec=5.00000 deg, Distance=1 kpc>'
+
+
+#<-------------------------Transformations------------------------------------->
+# Transformation functionality is the key to the whole scheme: they transform
+# low-level classes from one frame to another.
+
+# If no data (or `None`) is given, the class acts as a specifier of a frame, but
+# without any stored data.
+J2001 = astropy.time.Time('J2001',scale='utc')
+fk5_J2001_frame = coords.FK5(equinox=J2001)
+
+# if they do not have data, the string instead is the frame specification
+assert str(fk5_J2001_frame) == "<FK5 frame: equinox='J2000.000', obstime='B1950.000'>"
+
+#  Note that, although a frame object is immutable and can't have data added, it
+#  can be used to create a new object that does have data by giving the
+# `realize_frame` method a representation:
+srep = coords.SphericalRepresentation(lon=8*u.hour, lat=5*u.deg)
+fk5_j2001_with_data = fk5_J2001_frame.realize_frame(srep)
+assert fk5_j2001_with_data.data is not None
+# Now `fk5_j2001_with_data` is in the same frame as `fk5_J2001_frame`, but it
+# is an actual low-level coordinate, rather than a frame without data.
+
+# These frames are primarily useful for specifying what a coordinate should be
+# transformed *into*, as they are used by the `transform_to` method
+# E.g., this snippet precesses the point to the new equinox
+newfk5 = fk5.transform_to(fk5_J2001_frame)
+assert newfk5.equinox == J2001
+
+# classes can also be given to `transform_to`, which then uses the defaults for
+# the frame information:
+samefk5 = fk5.transform_to(coords.FK5)
+# `fk5` was initialized using default `obstime` and `equinox`, so:
+assert samefk5.ra == fk5.ra and samefk5.dec == fk5.dec
+
+
+# transforming to a new frame necessarily loses framespec information if that
+# information is not applicable to the new frame.  This means transforms are not
+# always round-trippable:
+fk5_2 =coords.FK5(ra=8*u.hour, dec=5*u.deg, equinox=J2001)
+ic_trans = fk5_2.transform_to(coords.ICRS)
+
+# `ic_trans` does not have an `equinox`, so now when we transform back to FK5,
+# it's a *different* RA and Dec
+fk5_trans = ic_trans.transform_to(coords.FK5)
+assert fk5_2.ra == fk5_trans.ra  # raises AssertionError
+
+# But if you explicitly give the right equinox, all is fine
+fk5_trans_2 = fk5_2.transform_to(coords.FK5(equinox=J2001))
+assert fk5_2.ra == fk5_trans_2.ra
+
+# Trying to tansforming a frame with no data is of course an error:
+coords.FK5(equinox=J2001).transform_to(coords.ICRS)  # raises ValueError
+
+
+# To actually define a new transformation, the same scheme as in the
+# 0.2/0.3 coordinates framework can be re-used - a graph of transform functions
+# connecting various coordinate classes together.  The main changes are:
+# 1) The transform functions now get the frame object they are transforming the
+#    current data into.
+# 2) Frames with additional information need to have a way to transform between
+#    objects of the same class, but with different framespecinfo values
+
+# An example transform function:
+@coords.dynamic_transform_matrix(SomeNewSystem, FK5)
+def new_to_fk5(newobj, fk5frame):
+    ot = newobj.obstime
+    eq = fk5frame.equinox
+    # ... build a *cartesian* transform matrix using `eq` that transforms from
+    # the `newobj` frame as observed at `ot` to FK5 an equinox `eq`
+    return matrix
+
+# Other options for transform functions include one that simply returns the new
+# coordinate object, and one that returns a cartesian matrix but does *not*
+# require `newobj` or `fk5frame` - this allows optimization of the transform.
+
+
+#<---------------------------"High-level" class-------------------------------->
+# The "high-level" class is intended to wrap the lower-level classes in such a
+# way that they can be round-tripped, as well as providing a variety of
+# convenience functionality.  This document is not intended to show *all* of the
+# possible high-level functionality, rather how the high-level classes are
+# initialized and interact with the low-level classes
+
+# this creates an object that contains an `ICRS` low-level class, initialized
+# identically to the first ICRS example further up.
+sc = coords.SkyCoordinate(coords.SphericalRepresentation(lon=8*u.hour, lat=5*u.deg, distance=1*u.kpc), system='icrs')
+# Other representations and `system` keywords delegate to the appropriate
+# low-level class. The already-existing registry for user-defined coordinates
+# will be used by `SkyCoordinate` to figure out what various the `system`
+# keyword actually means.
+
+# they can also be initialized using the preferred representation names
+sc = coords.SkyCoordinate(ra=8*u.hour, dec=5*u.deg, system='icrs')
+sc = coords.SkyCoordinate(l=120*u.deg, b=5*u.deg, system='galactic')
+
+# High-level classes can also be initialized directly from low-level objects
+sc = coords.SkyCoordinate(coords.ICRS(ra=8*u.hour, dec=5*u.deg))
+# The next example raises an error because the high-level class must always
+# have position data
+sc = coords.SkyCoordinate(coords.FK5(equinox=J2001))  # raises ValueError
+
+# similarly, the low-level object can always be accessed
+assert str(scoords.frame) == '<ICRS RA=120.000 deg, Dec=5.00000 deg>'
+
+# Should (eventually) support a variety of possible complex string formats
+sc = coords.SkyCoordinate('8h00m00s +5d00m00.0s', system='icrs')
+# In the next example, the unit is only needed b/c units are ambiguous.  In
+# general, we *never* accept ambiguity
+sc = coords.SkyCoordinate('8:00:00 +5:00:00.0', unit=(u.hour, u.deg), system='icrs')
+# The next one would yield length-2 array coordinates, because of the comma
+sc = coords.SkyCoordinate(['8h 5d', '2°5\'12.3" 0.3rad'], system='icrs')
+# It should also interpret common designation styles as a coordinate
+sc = coords.SkyCoordinate('SDSS J123456.89-012345.6', system='icrs')
+
+# the string representation can be inherited from the low-level class.
+assert str(sc) == '<SkyCoordinate (ICRS) RA=120.000 deg, Dec=5.00000 deg>'
+# but it should also be possible to provide formats for outputting to strings,
+# similar to `Time`.  This can be added right away or at a later date.
+
+# transformation is done the same as for low-level classes, which it delegates to
+scfk5_j2001 = scoords.transform_to(coords.FK5(equinox=J2001))
+
+# the key difference is that the high-level class remembers frame information
+# necessary for round-tripping, unlike the low-level classes:
+sc1 = coords.SkyCoordinate(ra=8*u.hour, dec=5*u.deg, equinox=J2001, system='fk5')
+sc2 = sc1.transform_to(coords.ICRS)
+# The next assertion succeeds, but it doesn't mean anything for ICRS, as ICRS
+# isn't defined in terms of an equinox
+assert sc2.equinox == J2001
+# But it *is* necessary once we transform to FK5
+sc3 = sc2.transform_to(coords.FK5)
+assert sc3.equinox == J2001
+assert sc1.ra == sc3.ra
+# Note that this did *not* work in the low-level class example shown above,
+# because the ICRS low-level class does not store `equinox`.
+
+
+# `SkyCoordinate` will also include the attribute-style access that is in the
+# v0.2/0.3 coordinate objects.  This will *not* be in the low-level classes
+sc = coords.SkyCoordinate(ra=8*u.hour, dec=5*u.deg, system='icrs')
+scgal = scoords.galactic
+assert str(scgal) == '<SkyCoordinate (Galactic) l=216.31707 deg, b=17.51990 deg>'
+
+
+# the existing `from_name` and `match_to_catalog_*` methods will be moved to the
+# high-level class as convenience functionality.
+
+m31icrs = SkyCoordinate.from_name('M31', system='icrs')
+assert str(m31icrs) == '<SkyCoordinate (ICRS) RA=10.68471 deg, Dec=41.26875 deg>'
+
+cat1 = SkyCoordinate(ra=[...]*u.hr, dec=[...]*u.deg, distance=[...]*u,kpc)
+cat2 = SkyCoordinate(ra=[...]*u.hr, dec=[...]*u.deg, distance=[...]*u,kpc
+idx2, sep2d, dist3d = cat1.match_to_catalog_sky(cat2)
+idx2, sep2d, dist3d = cat1.match_to_catalog_3d(cat2)
+
+
+# additional convenience functionality for the future should be added as methods
+# on `SkyCoordinate`, *not* the low-level classes.