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)
"""