Moon phase/illumination functions begun with PyEphem implementation - #39

A topic for discussion: looking at all of these functions that start with moon_*, do we really want to add these all as methods to the Observer class? I'm not sure either way, but it is starting to feel like the Observer class is getting a bit bloated...

@ejeschke Sorry, I missed it in the doctoring!

@bmorris3: I'm a bit confused by how this is now implemented. As I recall from one of your notebooks, the astropy/jplephem approach got the (GCRS) RA/Dec right (at least insofar as it matched pyephem), but it was only the GCRS->AltAz conversion that featured a bug. But the code you have here uses pyephem to compute the ra and dec but not the alt/az. Or am I misunderstanding something?

So I was thinking this would make it so that moon_altaz would still return a SkyCoord (in the AltAz frame), but where the alt/az was actually calulated using pyephem rather than the built-in astropy transformations.

@eteq You're right, the get_moon method produces the coordinates of the Moon without transforming to AltAz, since the GCRS->AltAz transformation is the buggy one. My intention is to mimic the get_sun function, which will return a SkyCoord with the RA/Dec of the Sun if computed with only a time as input, and if also computed with a location it will return the altaz. I've implemented the PyEphem and jplephem bits necessary to compute the RA/Dec without doing the AltAz transformation.

I (somehow) didn't realize you could initialize a SkyCoord from an AltAz frame, I guess I assumed it would be improperly defined without an associated location. I've updated the moon_altaz method to do so in d4aa918.

Ok gotcha. But one reason I was asking: calc_moon_phase_angle has an important subtlety: it's frame-dependent. That is, the angular separation of the moon from the sun is very different in a GCRS vs ICRS vs AltAz frame. In ICRS it's something very different from what you want because ICRS's origin is the solar system barycenter, rather than the Earth, so it's basically set by the earth and sun position relative to the barycenter, rather than the earth-moon-sun triangle. GCRS (which is what I think you're implicitly doing here) gets a lot closer, but it still is subtly different from what you'd see on the surface of the earth (which is what AltAz gives you).

Here, I think what you're doing in Observer.moon_phase is getting the moon as an AltAz, but the sun as GCRS. So I'd suggest adding one more step of converting the sun location to the appropriate AltAz before passing it into calc_moon_phase_angle.

@bmorris3 - I figured out what’s up with the recursion thing. get_moon currently returns the moon in the GCRS frame. So with the transform_to commented out, both sun and moon are in GCRS. So then sun.separation(moon) is just a direct geometric calculation, because the two are in the same frame. If you transform sun to AltAz, then moon will then has to get transformed to AltAz to do the comparison. But that can’t be done automatically because you haven’t given the moon SkyCoord a location. I’m not sure why that is causing an infinite recursion error rather than a helpful message about location being missing, but the root problem is the frames not matching.

So that brings in the reason for why moon_altaz being implemented via pyephem is important. The “best” way to solve this would instead be to give the moon location information so that it can just do the GCRS->AltAz transformation automatically. But that won’t work because of the bug that made us use pyephem.

So the work-around for now is to have get_moon yield something in the AltAz frame, rather than in GCRS. Then when you have the Sun in a matched AltAz frame, sun.separation(moon) will require no frame transformation and all is happy.

@eteq, I'm still getting the recursion depth error even after making the get_moon coordinate in the AltAz frame.

from astroplan.moon import get_moon
from astropy.coordinates import SkyCoord, EarthLocation, AltAz, get_sun
from astropy.time import Time
time = Time.now()
location = EarthLocation('-122d', '47d')
moon = get_moon(time, location)
sun = get_sun(time).transform_to(AltAz(location=location, obstime=time))
print(moon.separation(sun)) # Fails
print(sun.separation(moon)) # Fails

@bmorris3 - ok, now I think we've reached a bug in astropy rather than a problem with what frames are being compared. I suggest this as a work-around:
(this would replace the last two lines in the snippet you just showed)

from astropy.coordinates import Angle
from astropy.coordinates.angle_utilities import angular_separation

moonsph = moon.spherical
sunsph = moon.spherical

sep = Angle(angular_separation(moonsph.lon, moonsph.lat, sunsph.lon, sunsph.lat), u.deg)

that's basically the same as the guts of separation, but it short-circuits the frame transformation tests (which are the source of the problem)

@eteq / @cdeil I think this is as far as this PR is going to go until a few related astropy-related issues (astropy/astropy#4033, astropy/astropy#3998, astropy/astropy#3938, astropy/astropy#3920) are sorted out, and we can replace the PyEphem placeholders with the long-term jplephem equivalents.

I'm +1 to remove the calc_ in these functions:

• calc_moon_illumination -> moon_illumination
• calc_moon_phase_angle -> moon_phase_angle

Otherwise 👍

@cdeil I left the calc_s on because I was trying to keep calc_moon_illumination separate from the method on Observercalled moon_illumination. How can we best prevent overlapping names here?

I'm +0 on @cdeil's suggestion - it's fine to either remove the calc_ and otherwise leave it all as-is: from context the methods are always called as obs.moon_illumination while the function would be moon_illumination. And they do basically the same thing just with different information, right?

However, I think it's not that important either way, because those functions are not to be part of the public API, right? Rather, they're internal functions used by Observer.

Which reminds me of one additional (easy) change you need, @bmorris3: get_moon needs to be added to __all__ so that one can just do astroplan.get_moon instead of having to import the moon submodule.

If the function and method basically do the same thing, it's good if they have the same name, no?

Alright, looks good. Merging!