Great! I think it would be worthwhile implementing the correct 'tilt' between the Galactic coordinates frame and the galactic plane. The math is a pain, but @ChrisBeaumont, you have an implementation of that you can add here, right?

Actually I don't think it's that bad, right? Isn't it just a rotation about the y axis by

Or maybe I'm just being slow...

I need to check, it depends on the longitude you are looking at - if you look at l=0 then it's simple, but if l=90 then your line of sight is parallel to the plan, and if l=180 it looks above the plane.

@ChrisBeaumont - we had a repository with the transformation code somewhere right?

@adrn - I think you're right though, I think once you convert to cartesian then it's just an additional rotation around the y axis. Need to think about it.

@adrn - I think the reason we got confused before is because the issue is that b=0 doesn't cross the galactic center:

but crosses the plane behind the center! (plot copied from Alyssa's Nessie paper)

Oy....

But not to worry, @ChrisBeaumont has a magic function that does the conversion ;)

Guess it's basically a translation, then rotation, then translation back?

Thinking about this again a bit more, we should convert to cartesian coordinates in the frame of reference of the sun, then do a translation by the solar coordinates - it's just that the angle and rotation and translation will be such that b=0 doesn't line up with the GC.

Yeah, this is a huge pain (at least for my 3D-challenged brain). I eventually got some trustworthy IDL code from Erik Rosolowsky and Bob Benjamin, that I'm pretty sure does the math correctly. If nothing else, we could use this to add some unit tests.

I don't think there's any real standard definition of galactocentric coordinates, so you need to specify the height of the sun off the midplane and (to a lesser extent) the offset of (l,b,r)=(0,0,r_sun_gc) from SgrA* (nonzero in galactic coordinates)

So these numbers are slightly inconsistent with what's listed on wikipedia, but here are some example outputs from Bob Benjamin's code, assuming the sun is at (x, y, z) = (-8500, 0, 25) pc

I'm actually not sure that Sgr A*'s offset from the (l,b,d)=(0,0,Rsun) is accounted for in this

Galactocentric (x, y, z) points
           0             0             0
        1000             0             0
       -8500             0            30
       -8500             0            20
       -8600             0             0
       -9000         -1000            29

Galactic (l, b, d) points

       0.0000000000000000        -0.1685165033979717      8500.0367646263748611
       0.0000000000000000        -0.1507780108053125      9500.0328946798908873
       0.0000000000000000        89.9999950584399642         5.0000000000000000
       0.0000000000000000       -89.9999950584399642         5.0000000000000000
     179.9999901168799283       -14.0362426972493708       103.0776406404415155
    -116.5650447769313445         0.2049867266420032      1118.0411441445255605

You might also want to look at the Bessel maser parallax work, specifically

http://adsabs.harvard.edu/abs/2009ApJ...700..137R
and
http://adsabs.harvard.edu/abs/2014ApJ...783..130R

To answer @adrn's questions (IMHO):

Do we want to allow a transformation to a left-handed system as well? (because some people are crazy)

Boo crazy people! (i.e., lets leave it right-hand only for now unless someone requests left-hand)

Can we think of better names for the frame attributes?

I think they're great just the way they are.

One other thing, related to @ChrisBeaumont's comment about the Sgr A* offset: I don't think it necessarily makes sense to define this relative to the Galactic frame. The problem is that Galactic is just not a terribly well-defined frame (it was defined in 1960, after all!) and there's some intrinsic ambiguities about how to map it onto modern frames.

So I would instead suggest tying the system to ICRS. The simplest approach is to add three new frame attributes: galcen_ra, galcen_dec (in ICRS), and galcen_roll. Or, equivalently, galpole_ra, galpole_dec and lon0 (the latter is how the Galactic transform is done). Then just set the default to something that approximates Sgr A* along the x-axis and z-axis being the galactic pole.

Then to implement ICRS-> galactocentric, you just apply the three rotations before everything you have here (you can just copy-and-paste the code to create the matricies from the Galactic transform functions, sans the precession part). To go the other way, just invert it all.

Just saw this and wanted to say thanks for implementing this!
@adonath and I actually discussed whether it would be worth to implement this for gammapy/gammapy#157 yesterday.

This is failing with Numpy 1.5 on travis-ci:
https://travis-ci.org/astropy/astropy/jobs/30258541#L1532

But I think support for that will be dropped in the next stable release, so actually that travis-ci test should probably be dropped now so that new things in the Astropy development version like this don't have to support Numpy 1.5?

@eteq good idea -- added to my todo list, but there are some high priority things in the way right now...will get to it soon.

@adrn - it would be great to have this feature in Astropy - do you think you would have a chance to work on it soon? Maybe it doesn't matter exactly which convention is used at the moment as long as it's made clear. If you have a reference for the convention you go with, even better. Maybe in the long term we can even have different frames for different conventions?

I worked on this on a flight yesterday, so I'll try to finish it up this week

@adrn - any updates on this?

@eteq Oops, stuck on something minor that I need to chat with you about. You around to gitter or gchat tomorrow?

@adrn - sure, before 3:30 tomorrow would be fine. I'll log in to gchat and/or gitter when I arrive.

Without realizing this was here, I independently implemented the same thing:
https://gist.github.com/keflavich/83d047a15a6465c1ab52

@adrn @eteq - while chatting with @keflavich I came up with a diagram that summarizes different possible configurations - just so that we can always refer to configurations with a name - am I missing any?

I'm not saying which is 'right' at this point, but just trying to summarize the different point of views/implementations. I think this PR implements model 3, and I think @keflavich implements model 2.

@astrofrog - OK, yeah, I've been assuming 3 is the one we want. I'm not clear on how @keflavich is 2, though? Isn't his approach allowing for there to be an offset between (l,b)=(0,0) and his final x-axis?

Also, one question that @adrn and I got a bit stuck on : what are all the breakables necessary to provide here? It seems to me we need 5: to fix an arbitrary rotation from ICRS we need 3 variables, and then another 3 for the offset from the sun location to the GC. But we've also required that the GC-> sun axis be in the xz plane, so that removes 1 variable, leaving 5. Doors that sounds right?

I ask particularly now because @keflavich only has 4 variables (two angle for the GC direction in Galactic, zsun, and dsun). So what's the implicit extra constraint? (And do we want it?)

@astrofrog - OK, yeah, I've been assuming 3 is the one we want. I'm not clear on how @keflavich is 2, though? Isn't his approach allowing for there to be an offset between (l,b)=(0,0) and his final x-axis?

Strictly speaking, Model 2 is the correct one. Sgr A* is not a (l,b)=(0,0). In fact, b=0 intersects the true plane several kpc behind Sgr A* (12kpc from the sun), so it is not a small effect!

Also, one question that @adrn and I got a bit stuck on : what are all the breakables necessary to provide here? It seems to me we need 5: to fix an arbitrary rotation from ICRS we need 3 variables, and then another 3 for the offset from the sun location to the GC. But we've also required that the GC-> sun axis be in the xz plane, so that removes 1 variable, leaving 5. Doors that sounds right?

I think you only need 2 variables to specify an arbitrary rotation, right? I agree that we can assume the Sun-GC axis is in the xz plane (though note that this does not necessarily mean that the xz plane corresponds to l=0! (since Sgr A* is at l=359.9442)

I ask particularly now because @keflavich only has 4 variables (two angle for the GC direction in Galactic, zsun, and dsun). So what's the implicit extra constraint? (And do we want it?)

I'll double check this with @keflavich.

One question we need to answer - are we defining the Galactic center x=y=z=0 as Sgr A*? If so, then a good test of any transformation we come up with is that l=359.9442 b=-00.0462 should transform to (0,0,0).

@astrofrog @eteq on it! RE the Quantity issue -- should I just add the decompose() to make the warning go away?

weird! definitely a bug to report for Quantity, but decompose seems fine in this case.

Fun example to test out this PR - what is the latitude of points in the physical mid-plane?

import numpy as np
from astropy import units as u
from astropy.coordinates import Galactocentric, Galactic

x = np.linspace(-10., 10., 100) * u.kpc
y = np.linspace(-10., 10., 100) * u.kpc
X, Y = np.meshgrid(x, y)
Z = np.zeros(X.shape) * u.kpc

g = Galactocentric(x=X.ravel(), y=Y.ravel(), z=Z.ravel())
g2 = g.transform_to(Galactic)

lat = g2.b.degree.reshape(X.shape)

import matplotlib.pyplot as plt
fig = plt.figure()
ax = fig.add_subplot(1,1,1)
ax.imshow(lat, vmin=-1, vmax=+1, cmap=plt.cm.gist_heat, extent=[-10, 10, -10, 10], origin='lower')
cs = ax.contour(lat, levels=np.linspace(-0.5, 1, 16), colors='white', extent=[-10, 10, -10, 10])
ax.clabel(cs, cs.levels)
ax.set_xlabel('x (kpc)')
ax.set_ylabel('y (kpc)')
fig.savefig('latitude.png')

Similar shape to bottom panel (rotated 90 degrees) of:

https://www.authorea.com/users/23/articles/249/master/file/figures/nessie_faceon_both_final/nessie_faceon_both_final.png

This seems like it's doing the right thing, so I think that if you are also confident it's fine then we can merge it tonight and if we find any subtle bugs we still have time to fix them later (apart from the n-d issue reported above which we should fix now)

I've left two minor comments inline ... otherwise +1 to merge.

My example above now works without the ravel() - thanks!

👍 from me

@adrn - there is a real sphinx failure:

/home/travis/build/astropy/astropy/docs/api/astropy.coordinates.Galactocentric.rst:5: WARNING: py:obj reference target not found: l,b = (0º,0º)

/home/travis/build/astropy/astropy/docs/api/astropy.coordinates.Galactocentric.rst:5: WARNING: py:obj reference target not found: roll=0º

/home/travis/build/astropy/astropy/docs/api/astropy.coordinates.Galactocentric.rst:5: WARNING: py:obj reference target not found: l=90º

/home/travis/build/astropy/astropy/docs/api/astropy.coordinates.Galactocentric.rst:5: WARNING: py:obj reference target not found: b=90º

please use double backticks around these four bits. Since Travis is pretty busy today, it might be best to make sure it passes locally with python setup.py build_sphinx -w (this appears to hang for a while but it's because it's capturing stdout - just be patient).

@adrn - apart from the warnings this all looks good to me, so if you can fix these issues we can merge it and make the 1.0 deadline! :)

@astrofrog whoops -- fixed!

@adrn - just to be absolutely sure everything is fine, can you edit the last commit message and remove the [ci skip]? Thanks!

@astrofrog Oh, sorry, I thought you didn't want travis to run the build?

@adrn - I do in this case since the last build was a failure due to the doc errors. But otherwise this is good to go once we get green lights.

This can be merged once Travis and AppVeyor pass

Merged! Thanks @adrn