test-fix-gcrs.py

from __future__ import (absolute_import, division, print_function,
                        unicode_literals)

import numpy as np
import sys
sys.path.insert(0, '/Users/bmorris/git/tmp/astropy')
import astropy
from astropy.time import Time
from astropy.coordinates import (SkyCoord, Longitude, Latitude, GCRS,
                                 CartesianRepresentation, AltAz)
import astropy.units as u
from astropy.utils.data import download_file

def get_spk_file():
    """
    Get the Satellite Planet Kernel (SPK) file `de430.bsp` from NASA JPL.
    Download the file from the JPL webpage once and subsequently access a
    cached copy. This file is ~120 MB, and covers years ~1550-2650 CE [1]_.
    .. [1] http://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/planets/aareadme_de430-de431.txt
    """
    de430_url = ('http://naif.jpl.nasa.gov/pub/naif/'
                 'generic_kernels/spk/planets/de430.bsp')
    return download_file(de430_url, cache=True, show_progress=True)

def get_moon(time, location, pressure, use_pyephem=True):
    """
    SOURCE: https://github.com/astropy/astroplan/blob/edfed62c7034f655dcd40853d29432b14961182d/astroplan/moon.py
    
    Position of the Earth's moon.
    Currently uses PyEphem to calculate the position of the moon by default
    (requires PyEphem to be installed). Set ``use_pyephem`` to `False` to
    calculate the moon position with jplephem (requires jplephem to be
    installed).
    Parameters
    ----------
    time : `~astropy.time.Time` or see below
        Time of observation. This will be passed in as the first argument to
        the `~astropy.time.Time` initializer, so it can be anything that
        `~astropy.time.Time` will accept (including a `~astropy.time.Time`
        object).
    location : `~astropy.coordinates.EarthLocation`
        Location of the observer on Earth
    use_pyephem : bool (default = `True`)
        Calculate position of moon using PyEphem (requires PyEphem to be
        installed). If `False`, calculates position of moon with jplephem.
    Returns
    -------
    moon_sc : `~astropy.coordinates.SkyCoord`
        Position of the moon at ``time``
    """
    if not isinstance(time, Time):
        time = Time(time)

    if use_pyephem:
        try:
            import ephem
        except ImportError:
            raise ImportError("The get_moon function currently requires "
                              "PyEphem to compute the position of the moon.")

        moon = ephem.Moon()
        obs = ephem.Observer()
        obs.lat = location.latitude.to_string(u.deg, sep=':')
        obs.lon = location.longitude.to_string(u.deg, sep=':')
        obs.elevation = location.height.to(u.m).value
        if pressure is not None:
            obs.pressure = pressure.to(u.bar).value*1000.0

        if time.isscalar:
            obs.date = time.datetime
            moon.compute(obs)
            moon_alt = float(moon.alt)
            moon_az = float(moon.az)
            moon_distance = moon.earth_distance
        else:
            moon_alt = []
            moon_az = []
            moon_distance = []
            for t in time:
                obs.date = t.datetime
                moon.compute(obs)
                moon_alt.append(float(moon.alt))
                moon_az.append(float(moon.az))
                moon_distance.append(moon.earth_distance)
        return SkyCoord(alt=moon_alt*u.rad, az=moon_az*u.rad,
                        distance=moon_distance*u.AU,
                        frame=AltAz(location=location, obstime=time))
                
    else:
        try:
            import jplephem
        except ImportError:
            raise ImportError("The get_moon function currently requires "
                              "jplephem to compute the position of the moon.")

        # Calculate position of moon relative to Earth by subtracting the
        # vector pointing from the Earth-moon barycenter to the moon by
        # the vector from the Earth-moon barycenter to the Earth
        from jplephem.spk import SPK
        kernel = SPK.open(get_spk_file())
        cartesian_position = (kernel[3,301].compute(time.jd) -
                              kernel[3,399].compute(time.jd))
        x, y, z = cartesian_position*u.km

        # Convert to GCRS coordinates
        cartrep = CartesianRepresentation(x=x, y=y, z=z)
        return SkyCoord(cartrep, frame=GCRS(obstime=time,
                                            obsgeoloc=location))


from astroplan import get_site, Observer

time = Time('1984-01-21 12:00:00')
location = get_site('APO')
obs = Observer(location=location, pressure=0)

moon_pyephem = get_moon(time, location, 0*u.bar)
moon_jplephem = get_moon(time, location, 0*u.bar, use_pyephem=False)

altaz_pyephem = obs.altaz(time, moon_pyephem)
altaz_jplephem = obs.altaz(time, moon_jplephem)

print("With astropy v{0}".format(astropy.__version__))
print("PyEphem (Alt, Az): ({0.alt}, {0.az})".format(altaz_pyephem))
#print("\t (RA, Dec): ({0.ra}, {0.dec})".format(moon_pyephem))

print("\njplephem (Alt, Az): ({0.alt}, {0.az})".format(altaz_jplephem))
#print("\t (RA, Dec): ({0.ra}, {0.dec})".format(moon_jplephem))

"""
Output with astropy v1.0.5
PyEphem (Alt, Az): (55.4920153458 deg, 239.031895666 deg)
jplephem (Alt, Az): (34.177292986 deg, 135.374232837 deg)
"""