Conversion from celestial to altaz coordinates

[mperrin]

At the #hackAAS hack day, I started helping @elovegrove to modify her program PYRO (https://github.com/elovegrove/PYRO) to use astropy.coordinates. However I quickly discovered that astropy.coordinates is still missing an interface to convert from celestial coordinates to apparent horizontal coordinates at some observatory. So, I started coding.

The first ingredient is some way to represent an observer's position on the earth. I implemented a LatitudeLongitude coordinates class, following the examples of the other systems. For now there's only one such system and it doesn't have any transformations, and also it ignores the details of the geoid to just treat the Earth as a perfect sphere. (I imagine if eventually more precision is needed, one could use the transformations framework to define conversions between various geoid models… but that's getting ahead of the story).

Then, the next question is, what should the API be for using LatitudeLongitudes to get a HorizontalCoordinates? I propose something like this:

>>> kpno = coord.LatitudeLongitude(31.9583, -111.5967, unit=(u.degree, u.degree))
>>> star = coord.ICRSCoordinates(ra=10.68458, dec=41.26917, unit=(u.degree, u.degree))
>>> time = astropy.time.Time('2014-01-15 01:02:03', scale='utc')
>>>
>>> star.observed_from(kpno, time)
<HorizontalCoordinates object>

One could conversely imagine making it a function on the observatory location instead:

>>> kpno.observe(star, time)
<HorizontalCoordinates object>

Either API could work fine. It just seems a bit cleaner to me to have the time and place of the observation treated identically, as part of the function call. It should still support a syntax using the convert_to function, something like

>>> altaz = c.convert_to(HorizontalCoordinate, location=kpno, time=obstime)

but I think this is a common enough conversion that it deserves some syntactic sugar. Thoughts?

[timj]

Some of this came up in the initial API discussions for coordinates. In particular the need for time and observer location and whether those are top level attributes of all coordinate objects or something that is added to subclasses. The other discussion concerned whether NDData should have epoch information at the top level or rely on it being in the coordinates object.

The pyast library handles AZEL. @dsberry had implemented a version of the astropy coordinates API using pyast.

I think the next step is to define the observer object. I'm guessing that spacecraft as observers causes issues.

[astrofrog]

Just a quick comment, I think that the position on the earth should just be represented by a Location class as there are other things that one could add, e.g. altitude.

[astrofrog]

A follow-up on my previous comment - one could simply have Location as a base class for any 3-d observer position, which could be a spacecraft, observatory, etc. Then, EarthLocation could be defined with a longitude, latitude, and optionally altitude. I'm not saying all should be implemented straight away, but I think such an API would make sense.

[cdeil]

This is such a complex topic, because there's so many different ways to make this API, i.e. how to group the different parameters needed into different classes.

Here's a (depending on the scope of astropy.coordinates incomplete) list of parameters:

• coordinate lon and lat and epoch
• time of observation
• observer lon and lat and elevation
• observer temp and pressure and wavelength if we want astropy to handle refraction

At the moment we have Time and Coordinate objects already and it's clear that we need another object, which could be a Location or an Observer, which is what pyephem has which as far as I know is the only existing easy-to-use Python package for celestial <=> horizontal coordinate conversions.

One issue I think we should keep an eye on is that we don't duplicate the same parameters in different objects, or at least as few duplications as possible if it can't be avoided. E.g. currently astropy.time.Time has observer lon and lat as attributes! And astropy.coordinates.FK5Coordinates has has equinox and obstime as Time attributes, each with their own equinox.lon, equinox.lat, obstime.lon, obstime.lat attributes. Now in the API you propose above you have the obstime again as a standalone object, i.e. when you call

kpno.observe(star, time)

and star.obstime and time are different, do you throw an error or do you ignore star.obstime?

Another point I'd like to make is that I think it would be a mistake to implement celestial <=> horizontal coordinate conversions from scratch, because it's a huge effort. Instead we should wrap an existing high-precision and fast C library and give it a nice Python API (does anyone disagree on this point?).

A while back I've tried to summarize the existing libraries and tried to identify which would have a BSD-compatible license here. Unfortunately as far as I can tell there is no BSD-compatible suitable library available at the moment, so we should bug the authors of the candidate libraries to re-consider:

• SLALIB and AST are GPL-licensed (could ask them to re-consider?)
• What is the TPM license?
• An the astropy-dev mailing list there has been discussion about the SOFA library, as far as I know we did ask them to re-licese as BSD-compatible, but there was no reply.

To get a taste for what is needed for high-precision celestial <=> horizontal coordinate conversions, have a look at the TPM states and the TPM state parameters. I don't see the value of implementing a very low-precision conversion assuming the Earth is a sphere as suggested above (anyone disagree?)

I'm not sure, but I doubt that the API we are discussing here is independent of the library we base it on, so I doubt it is a good idea to implement some nice API for low-precision coordinate conversions for astropy 0.3 and assume this API would work well (clean and efficient) for later when we do change to fast, high-precision coordinate conversions using some library. Which parameters should be stored in which astropy objects depends on the scope and precision we eventually want to achieve, e.g. TPM uses 6-vector coordinates (3D cartesian position and velocity) as the basic object, but as far as I know no-one is considering velocity vectors for astropy.coordinates at the moment, and I simply don't know if the precision without them will be good enough for astropy or if it would be awkward / partly incorrect to use TPM by always dropping the velocity part of the 6-vectors.

So before or in parallel to this discussion on the astropy.coordinates API, I propose we start discussing these questions:

• Which C libraries (TPM, AST, SLALIB, SOFA, ...?) are options to base astropy.coordinates on?
• What precision / scope do we want for astropy.coordinates?
• In the astropy.coordinates API, which attributes can be arrays and which have to be scalars, i.e. which use-cases will be fast (both because the underlying library doesn't recompute everything (like nutation, ...) for each array element or because the loop over array elements occurs in C land or Python land)?

For the last two questions, let me list some of the coordinate computations I'd like to be able to do (with high precision and if possible efficiently) with astropy.coordinates:

• Compute the (Alt, Az) of my telescope has to track to follow a given star in a 1 hour observation: the coordinate and the observer position are scalars and the observation time is an array.
• Compute the (Alt, Az) for a given list of a few 1000 observations spread over a few years, i.e. compared to the previous case the equinox is now different for each obstime (e.g. TPM has extra states for this transformation from Heliocentric mean celestial coordinates to Topocentric apparent coordinates) and possibly also the pressure and temperature for each observation is known and different.
• Compute the (Alt, Az) for a catalog of stars, the observer and time are scalars, and the coordinate is an array.
• Compute the ICRS (RA, DEC) for arrays of (ALT, AZ, Time) and a fixed observer
• Compute the (Alt, Az, Distance) for arrays of points in the Earth's atmosphere, i.e. given arrays of (X, Y, Z, Time) and a fixed observer
• Atmospheric refraction effects

These things I would say are out of the scope for astropy.coordinates (anyone think they should go in the astropy core?):

• ephemerides (sun, moon, planets, Earth orbit, ...)
• outer space observers (in a spacecraft, on Moon, Mars, Vega, M31)
  *Barycentric corrections

Sorry for the long post!

I guess my comment can be summarized like this:
Yes, we need horizontal coordinates and a Location and maybe also Observer class for astropy.coordinates, but it would make sense to discuss the "big picture" or "vision" or "scope" for astropy.coordinates first before discussing the API.

[eteq]

@mperrin - To start with a minor detail: I think it might be better to start this as a new document (observed_coordinates.py or something) instead of modifying the existing coordinates.py - the intent was for these documents to be left mostly as-is and new significant APIs get their own files.

@mperrin @astrofrog - I definitely like what you're saying here, and is exactly what I was thinking. We can combine the two by having a Location include a LatitudeLongitudes inside it.

I also think that adding this functionality should also add a registry of sites, starting with some of the "standard" lists of observatories. Later on, we would expand this to include additional information (or at least the option of it) like telescopes, various parameters of those telescopes, instruments on them, and even information like particular CCDs' readout characteristics. We don't want to be in the business of managing all this information in the source code, but it would be relatively easy to make a framework to do this so that astropy data reduction/analysis tools can make use of whatever a user supplies.

As for @mperrin's question of ICRSCoordinates.observed_from() vs. Location.observe(), I have no strong opinion. My inclination is actually to do both, although only implement it in one of them (probably the former) and call it from the other. I think there are a variety of cases where either are more convenient than the other. But @cdeil is definitely right that this requires care to make sure information isn't duplicated across the two.

Speaking of which, the particular case of kpno.observe(star, time) is probably the most important one. My inclination is that it should be kpno.observe(star), and if there's no obstime in the star, it raises an exception that says you need to give it one. There could be optional syntactic sugar kpno.observe(star, time) that temporarily sets obstime, but it should issue a warning if there's already an obstime there (or maybe even an exception?). All the other cases of "additional information" should be dealt with the same way: pick which one it should live in, and the API only says you access it from there, but other functions/methods may temporarily fiddle with it (although never silently overriding something already there).

[eteq]

@timj @astrofrog re: space-based locations: Actually, I would think these are relatively easy. The complications for earth coordinates mostly come from dealing with equinoxes and such, which are not an issue for spacecraft observers. If you have e.g. barycentric coordinates (cartesian or otherwise) for your spacecraft, you can take any ICRSCoordinate, offset the origin to the location of your spacecraft, possibly rotate, switch back to spherical, and you're done. That machinery is already in astropy.coordinates, and in fact most of the necessary steps are already coded up in a similar way in astropysics.

(This post is not GR-compliant, but I'm willing to punt that for later... :)

[eteq]

@cdeil - A discussion of what's "in scope" is exactly what the API discussion is about, isn't it? I do not want this to be in general terms, though, as that's what many of these discussions have devolved into. So I like that @mperrin has a specific set of requirements/examples here, and I think I agree with all your items of what should be in (although I think atmospheric refraction is farther down on the todo list than all the others).

I think ephemerides should be in Astropy - not a huge set of them, but a framework that others can build from if they want to. We definitely need Sun and Moon (and all the hard work is already done, as that's in SOFA, which we already have), and I think the planets wouldn't be too much trouble. It also wouldn't be that much work to add in a way to query the JPL ephemerides files and represent them in a way Astropy could read. I think this should all be in the core, but is the next item on the list after getting the coordinates hashed out. (Also, as I mentioned above, outer space observers should actually be pretty easy to implement in the existing framework.) We should also keep possible future compatibility with GR in mind, but it's last on the priority list (but yes, in the core someday is a future goal).

I doubt that the API we are discussing here is independent of the library we base it on

I'm going to push back a little here: I strongly disagree that we should be thinking about the C library implementations when designing the API. We should design an API that makes sense, and then implement accordingly. I think a lot of the other library's APIs are confusing because they did not do this, or added on API features once they realized they could implement them.

I have plenty I could say about the rest, but I don't want to get caught up here in the question of whether or not to use a C library, as it's a separate discussion. There will definitely be a skype/google hangour-con on this topic soon after 0.2 is released - I presume you want to be in on that @cdeil ?

[mperrin]

@eteq - re "My inclination is that it should be kpno.observe(star), and if there's no obstime in the star"

The whole notion of having an "obstime in the star" sounds painfully unintuitive to me. That's not a meaningful concept physically, and it's not how most astronomers think about this topic. Yes, there are some coordinates types that include an obstime parameter, but that's really the exception not the rule.

Consider the most typical case: You've got a star which has some known ICRS/J2000 coordinates, and you want to figure out its position on some other datetime, say today 2013 Jan 21 10:00 pm EST. Thus the datetime of the observation is completely different from the datetime associated with the coordinates class. Trying to conflate those is going to cause far more confusion than anything else, in my opinion.

[mperrin]

Re: space based observers, this is not quite as simple as you make it out to be. The relative-motion-induced shift in apparent locations ("aberration of starlight") is a larger effect for space based observers than terrestrial ones given the higher velocities. You can't actually punt this effect without getting relative positions wrong by many tens of arcseconds.

That said, I think we should (for now, and likely a long time) set aside space based observers entirely. The fundamental paradigm of an observer with a fixed location and a locally-defined horizontal coordinate system doesn't apply. Spacecraft attitude control is a complex and subtle subject which I don't see astropy making inroads into any time soon. I'd suggest waiting until we have some clear use cases for this before giving it any more thought. We have plenty of well-defined and relevant use cases for terrestrial calculations to keep us busy for now. :-)

[timj]

@eteq The list of telescope coordinates is a pretty easy thing to deal with. The perl Astro::Telescope (timj/perl-Astro-Telescope@01abf3d) module handles it by supporting the names of pre-canned coordinates that are available in SLALIB or PAL but it also generates coordinates from the list at the Minor Planet Center (http://www.minorplanetcenter.net/iau/lists/ObsCodes.html).

[astrofrog]

@eteq - we may want to consider whether the obstime in the Coordinates object is really the same as the observation time a user wants - obstime in the coordinates is really the date/time when the coordinates were measured and are valid (for moving frames), which might be different from the observation time the user wants now.

For example, one could use FK4 coordinates with an equinox of B1950, and epoch (for the coordinates) of B1978, and want to observe that source tomorrow. The user can't just update the obstime attribute, otherwise the coordinates are wrong. They need to transform the coordinates to the new observation date/time first before transforming to alt/az, and this is not really obvious.

So one could rename obstime to something else, and in most cases would not be needed (only for frames like FK4 where it makes a big difference), and then the date/time of the desired observation could be passed into the observe() function or method. Otherwise, if a user wants to find Alt/Az vs time for a given coordinate object, they will have to update the obstime attribute many times, which doesn't seem like the right way to do things. If we decouple the date/time of definition of the coordinates (currently obstime) and the date/time of observation, which is not an attribute, this is no longer an issue.

[timj]

@eteq regarding plannets and moon ephemeris data the SOFA implementations are fine for low precision situations (up to an arcminute error I think). There are a number of other approaches that are more accurate but without requiring the full DE405 JPL ephemeris data files. Pat Wallace suggested that we take a look at http://cdsarc.u-strasbg.fr/viz-bin/Cat?VI/87 since it gives milliarcsec accuracy and only requires 1 MB of data file to be included. I've also just discovered Chapront 1995, A&AS, 109 which is used by some telescope control systems but I haven't had a look at the size of the data files required to implement it.

[peterjc]

For instrumenting a telescope I was just looking for sample code using astropy for conversion between celestial coordinates (RA and declination) and local Earth based observer's horizon based Altitude/Azimuth coordinates (taking into account the observer's location - latitude, longitude, and perhaps altitude, and the current date/time). [In both directions ideally]

I didn't find anything obvious in the documentation, and this pull request suggests it is an open feature request... is that correct? Thanks!

[astrofrog]

@peterjc - indeed, it is not implemented yet, as there is some initial refactoring of astropy.coordinates scheduled first (cc @eteq)