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.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))
    elif use_skyfield:
        from skyfield import api

        if use_skyfield is True:
            use_skyfield = 'de421.bsp'
        planets = api.load(use_skyfield)
        loc = planets['earth'].topos(latitude_degrees=location.latitude.deg,
                                     longitude_degrees=location.longitude.deg)
        moon_obs = loc.at(api.JulianDate(tai=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:
        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())
        if use_icrs:
            cartesian_position = (kernel[0,3].compute(time.jd) +
                                  kernel[3,301].compute(time.jd))
            cartrep = CartesianRepresentation(cartesian_position*u.km)
            moon_icrs = SkyCoord(cartrep, frame='icrs')
            return moon_icrs.transform_to(AltAz(obstime=time, location=location))
        else:
            cartesian_position = (kernel[3,301].compute(time.jd) -
                                  kernel[3,399].compute(time.jd))

            # Convert to GCRS coordinates
            cartrep = CartesianRepresentation(cartesian_position*u.km)
            # return SkyCoord(cartrep, frame=GCRS(obstime=time,
            #                                     obsgeoloc=location))
            return SkyCoord(cartrep, frame=GCRS(obstime=time,
                                                obsgeoloc=u.Quantity(location.to_geocentric())))


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, 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)
moon_jplephem_icrs = get_moon(time, location, 0*u.bar, use_icrs=True)

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

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))
print("jplephem/Astropy(ICRS) (Alt, Az, Dist): ({0.alt}, {0.az}, {0.distance.km} km)".format(altaz_jplephem_icrs))

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

Trying to execute your code in Spyder (Anaconda).

[#################################] 100% de421.bsp
Traceback (most recent call last):

File "C:\Anaconda2\lib\site-packages\IPython\core\interactiveshell.py", line 3083, in run_code
self.showtraceback()

File "C:\Anaconda2\lib\site-packages\IPython\core\interactiveshell.py", line 1880, in showtraceback
value, tb, tb_offset=tb_offset)

File "C:\Anaconda2\lib\site-packages\IPython\core\ultratb.py", line 1242, in structured_traceback
self, etype, value, tb, tb_offset, number_of_lines_of_context)

File "C:\Anaconda2\lib\site-packages\IPython\core\ultratb.py", line 1159, in structured_traceback
self, etype, value, elist, tb_offset, number_of_lines_of_context

File "C:\Anaconda2\lib\site-packages\IPython\core\ultratb.py", line 511, in structured_traceback
lines = ''.join(self._format_exception_only(etype, value))

File "C:\Anaconda2\lib\site-packages\IPython\core\ultratb.py", line 623, in _format_exception_only
Colors.Normal, s))

UnicodeDecodeError: 'ascii' codec can't decode byte 0xcf in position 11: ordinal not in range(128)

I also meet the same error