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, ICRS,
                                 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=False, use_skyfield=False, use_icrs=False):
    """
    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.lat.to_string(u.deg, sep=':')
        obs.lon = location.lon.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))
    elif use_skyfield:
        from skyfield import api

        if use_skyfield is True:
            use_skyfield = 'de421.bsp'
        planets = api.load(use_skyfield)
        ts = api.load.timescale()
        loc = planets['earth'] + api.Topos(latitude_degrees=location.lat.deg,
                                     longitude_degrees=location.lon.deg)
        moon_obs = loc.at(ts.tai(jd=time.tai.jd)).observe(planets['moon'])
        alt, az, d = moon_obs.apparent().altaz()
        return SkyCoord(alt=alt.radians*u.rad, az=az.radians*u.rad,
                        distance=d.au*u.AU,
                        frame=AltAz(location=location, obstime=time))



    else:
        return astropy.coordinates.get_body('moon', time, location, 'de430').transform_to(AltAz(obstime=time, location=location))


from astroplan import Observer
from astropy.coordinates import EarthLocation

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

moon_pyephem = get_moon(time, location, 0*u.bar, use_pyephem=True)
moon_skyfield = get_moon(time, location, 0*u.bar, use_skyfield=True)
moon_skyfield430 = get_moon(time, location, 0*u.bar, use_skyfield='de430.bsp')
moon_jplephem = get_moon(time, location, 0*u.bar)

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, Dist): ({0.alt}, {0.az}, {0.distance.km} km)".format(altaz_pyephem))

print("Skyfield (Alt, Az, Dist): ({0.alt}, {0.az}, {0.distance.km} km)".format(moon_skyfield))
print("Skyfield(de430) (Alt, Az, Dist): ({0.alt}, {0.az}, {0.distance.km} km)".format(moon_skyfield430))

print("jplephem/Astropy (Alt, Az, Dist): ({0.alt}, {0.az}, {0.distance.km} km)".format(altaz_jplephem))

$ pip list | grep -E 'astropy|astroplan|skyfield|pyephem'
astroplan                0.7
astropy                  4.3.dev234+g2c6fd69b4 /home/juanlu/Personal/Astropy/astropy
pyephem                  3.7.7.0
pytest-astropy           0.8.0
pytest-astropy-header    0.1.2
skyfield                 1.33
$ python test_fix_gcrs.py 
With astropy v4.3.dev234+g2c6fd69b4
PyEphem (Alt, Az, Dist): (55.492015345752876 deg, 239.0318956655817 deg, 356424.8481068663 km)
Skyfield (Alt, Az, Dist): (55.49124334373045 deg, 239.03351681046558 deg, 356403.6954445954 km)
Skyfield(de430) (Alt, Az, Dist): (55.49124352738164 deg, 239.03351684329135 deg, 356403.69525128463 km)
jplephem/Astropy (Alt, Az, Dist): (56.48500486073724 deg, 238.7107927029978 deg, 357213.21025216894 km)
jplephem/Astropy(ICRS) (Alt, Az, Dist): (55.41730784517011 deg, 239.37137207735177 deg, 355352.7959745067 km)

With latest Skyfield and corrections:

$ python test_fix_gcrs.py 
With astropy v4.3.dev234+g2c6fd69b4
PyEphem (Alt, Az, Dist): (55.492015345752876 deg, 239.0318956655817 deg, 356424.8481068663 km)
Skyfield (Alt, Az, Dist): (55.49124334373045 deg, 239.03351681046558 deg, 356403.6954445954 km)
Skyfield(de430) (Alt, Az, Dist): (55.49124352738164 deg, 239.03351684329135 deg, 356403.69525128463 km)
jplephem/Astropy (Alt, Az, Dist): (55.49118251346562 deg, 239.03406534146313 deg, 356399.58247342776 km)