# Licensed under a 3-clause BSD style license - see LICENSE.rst
from __future__ import (absolute_import, division, print_function,
unicode_literals)
from six import string_types
# Standard library
import sys
import datetime
import warnings
# Third-party
from astropy.coordinates import (EarthLocation, SkyCoord, AltAz, get_sun,
get_moon, Angle, Longitude)
import astropy.units as u
from astropy.time import Time
from astropy.utils.exceptions import AstropyDeprecationWarning
import numpy as np
import pytz
# Package
from .exceptions import TargetNeverUpWarning, TargetAlwaysUpWarning
from .moon import moon_illumination, moon_phase_angle
from .target import get_skycoord, SunFlag, MoonFlag
__all__ = ["Observer"]
MAGIC_TIME = Time(-999, format='jd')
# Handle deprecated MAGIC_TIME variable
def deprecation_wrap_module(mod, deprecated):
"""Return a wrapped object that warns about deprecated accesses"""
deprecated = set(deprecated)
class DeprecateWrapper(object):
def __getattr__(self, attr):
if attr in deprecated:
warnmsg = ("`MAGIC_TIME` will be deprecated in future versions "
"of astroplan. Use masked Time objects instead.")
warnings.warn(warnmsg, AstropyDeprecationWarning)
return getattr(mod, attr)
return DeprecateWrapper()
sys.modules[__name__] = deprecation_wrap_module(sys.modules[__name__],
deprecated=['MAGIC_TIME'])
def _generate_24hr_grid(t0, start, end, n_grid_points, for_deriv=False):
"""
Generate a nearly linearly spaced grid of time durations.
The midpoints of these grid points will span times from ``t0``+``start``
to ``t0``+``end``, including the end points, which is useful when taking
numerical derivatives.
Parameters
----------
t0 : `~astropy.time.Time`
Time queried for, grid will be built from or up to this time.
start : float
Number of days before/after ``t0`` to start the grid.
end : float
Number of days before/after ``t0`` to end the grid.
n_grid_points : int (optional)
Number of grid points to generate
for_deriv : bool
Generate time series for taking numerical derivative (modify
bounds)?
Returns
-------
`~astropy.time.Time`
"""
if for_deriv:
time_grid = np.concatenate([[start - 1 / (n_grid_points - 1)],
np.linspace(start, end,
n_grid_points)[1:-1],
[end + 1 / (n_grid_points - 1)]]) * u.day
else:
time_grid = np.linspace(start, end, n_grid_points) * u.day
# broadcast so grid is first index, and remaining shape of t0
# falls in later indices. e.g. if t0 is shape (10), time_grid
# will be shape (N, 10). If t0 is shape (5, 2), time_grid is (N, 5, 2)
while time_grid.ndim <= t0.ndim:
time_grid = time_grid[:, np.newaxis]
# we want to avoid 1D grids since we always want to broadcast against targets
if time_grid.ndim == 1:
time_grid = time_grid[:, np.newaxis]
return t0 + time_grid
[docs]class Observer(object):
"""
A container class for information about an observer's location and
environment.
Examples
--------
We can create an observer at Subaru Observatory in Hawaii two ways. First,
locations for some observatories are stored in astroplan, and these can be
accessed by name, like so:
>>> from astroplan import Observer
>>> subaru = Observer.at_site("Subaru", timezone="US/Hawaii")
To find out which observatories can be accessed by name, check out
`~astropy.coordinates.EarthLocation.get_site_names`.
Next, you can initialize an observer by specifying the location with
`~astropy.coordinates.EarthLocation`:
>>> from astropy.coordinates import EarthLocation
>>> import astropy.units as u
>>> location = EarthLocation.from_geodetic(-155.4761*u.deg, 19.825*u.deg,
... 4139*u.m)
>>> subaru = Observer(location=location, name="Subaru", timezone="US/Hawaii")
You can also create an observer without an
`~astropy.coordinates.EarthLocation`:
>>> from astroplan import Observer
>>> import astropy.units as u
>>> subaru = Observer(longitude=-155.4761*u.deg, latitude=19.825*u.deg,
... elevation=0*u.m, name="Subaru", timezone="US/Hawaii")
"""
@u.quantity_input(elevation=u.m)
def __init__(self, location=None, timezone='UTC', name=None, latitude=None,
longitude=None, elevation=0*u.m, pressure=None,
relative_humidity=None, temperature=None, description=None):
"""
Parameters
----------
location : `~astropy.coordinates.EarthLocation`
The location (latitude, longitude, elevation) of the observatory.
timezone : str or `datetime.tzinfo` (optional)
The local timezone to assume. If a string, it will be passed
through ``pytz.timezone()`` to produce the timezone object.
name : str
A short name for the telescope, observatory or location.
latitude : float, str, `~astropy.units.Quantity` (optional)
The latitude of the observing location. Should be valid input for
initializing a `~astropy.coordinates.Latitude` object.
longitude : float, str, `~astropy.units.Quantity` (optional)
The longitude of the observing location. Should be valid input for
initializing a `~astropy.coordinates.Longitude` object.
elevation : `~astropy.units.Quantity` (optional), default = 0 meters
The elevation of the observing location, with respect to sea
level. Defaults to sea level.
pressure : `~astropy.units.Quantity` (optional)
The ambient pressure. Defaults to zero (i.e. no atmosphere).
relative_humidity : float (optional)
The ambient relative humidity.
temperature : `~astropy.units.Quantity` (optional)
The ambient temperature.
description : str (optional)
A short description of the telescope, observatory or observing
location.
"""
self.name = name
self.pressure = pressure
self.temperature = temperature
self.relative_humidity = relative_humidity
# If lat/long given instead of EarthLocation, convert them
# to EarthLocation
if location is None and (latitude is not None and
longitude is not None):
self.location = EarthLocation.from_geodetic(longitude, latitude,
elevation)
elif isinstance(location, EarthLocation):
self.location = location
else:
raise TypeError('Observatory location must be specified with '
'either (1) an instance of '
'astropy.coordinates.EarthLocation or (2) '
'latitude and longitude in degrees as '
'accepted by astropy.coordinates.Latitude and '
'astropy.coordinates.Latitude.')
# Accept various timezone inputs, default to UTC
if isinstance(timezone, datetime.tzinfo):
self.timezone = timezone
elif isinstance(timezone, string_types):
self.timezone = pytz.timezone(timezone)
else:
raise TypeError('timezone keyword should be a string, or an '
'instance of datetime.tzinfo')
def __repr__(self):
"""
String representation of the `~astroplan.Observer` object.
Examples
--------
>>> from astroplan import Observer
>>> keck = Observer.at_site("Keck", timezone="US/Hawaii")
>>> print(keck) # doctest: +FLOAT_CMP
<Observer: name='Keck',
location (lon, lat, el)=(-155.478333333 deg, 19.8283333333 deg, 4160.0 m),
timezone=<DstTzInfo 'US/Hawaii' LMT-1 day, 13:29:00 STD>>
"""
class_name = self.__class__.__name__
attr_names = ['name', 'location', 'timezone', 'pressure', 'temperature',
'relative_humidity']
attr_values = [getattr(self, attr) for attr in attr_names]
attributes_strings = []
for name, value in zip(attr_names, attr_values):
if value is not None:
# Format location for easy readability
if name == 'location':
formatted_loc = ["{} {}".format(i.value, i.unit)
for i in value.to_geodetic()]
attributes_strings.append(
"{} (lon, lat, el)=({})".format(
name, ", ".join(formatted_loc)))
else:
if name != 'name':
value = repr(value)
else:
value = "'{}'".format(value)
attributes_strings.append("{}={}".format(name, value))
return "<{}: {}>".format(class_name, ",\n ".join(attributes_strings))
[docs] @classmethod
def at_site(cls, site_name, **kwargs):
"""
Initialize an `~astroplan.observer.Observer` object with a site name.
Extra keyword arguments are passed to the `~astroplan.Observer`
constructor (see `~astroplan.Observer` for available keyword
arguments).
Parameters
----------
site_name : str
Observatory name, must be resolvable with
`~astropy.coordinates.EarthLocation.get_site_names`.
Returns
-------
`~astroplan.observer.Observer`
Observer object.
Examples
--------
Initialize an observer at Kitt Peak National Observatory:
>>> from astroplan import Observer
>>> import astropy.units as u
>>> kpno_generic = Observer.at_site('kpno')
>>> kpno_today = Observer.at_site('kpno', pressure=1*u.bar, temperature=0*u.deg_C)
"""
name = kwargs.pop('name', site_name)
if 'location' in kwargs:
raise ValueError("Location kwarg should not be used if "
"initializing an Observer with Observer.at_site()")
return cls(location=EarthLocation.of_site(site_name), name=name, **kwargs)
[docs] def astropy_time_to_datetime(self, astropy_time):
"""
Convert the `~astropy.time.Time` object ``astropy_time`` to a
localized `~datetime.datetime` object.
Timezones localized with `pytz`_.
.. _pytz: https://pypi.python.org/pypi/pytz/
Parameters
----------
astropy_time : `~astropy.time.Time`
Scalar or list-like.
Returns
-------
`~datetime.datetime`
Localized datetime, where the timezone of the datetime is
set by the ``timezone`` keyword argument of the
`~astroplan.Observer` constructor.
Examples
--------
Convert an astropy time to a localized `~datetime.datetime`:
>>> from astroplan import Observer
>>> from astropy.time import Time
>>> subaru = Observer.at_site("Subaru", timezone="US/Hawaii")
>>> astropy_time = Time('1999-12-31 06:00:00')
>>> print(subaru.astropy_time_to_datetime(astropy_time))
1999-12-30 20:00:00-10:00
"""
if not astropy_time.isscalar:
return [self.astropy_time_to_datetime(t) for t in astropy_time]
# Convert astropy.time.Time to a UTC localized datetime (aware)
utc_datetime = pytz.utc.localize(astropy_time.utc.datetime)
# Convert UTC to local timezone
return self.timezone.normalize(utc_datetime)
[docs] def datetime_to_astropy_time(self, date_time):
"""
Convert the `~datetime.datetime` object ``date_time`` to a
`~astropy.time.Time` object.
Timezones localized with `pytz`_. If the ``date_time`` is naive, the
implied timezone is the ``timezone`` structure of ``self``.
Parameters
----------
date_time : `~datetime.datetime` or list-like
Returns
-------
`~astropy.time.Time`
Astropy time object (no timezone information preserved).
Examples
--------
Convert a localized `~datetime.datetime` to a `~astropy.time.Time`
object. Non-localized datetimes are assumed to be UTC.
<Time object: scale='utc' format='datetime' value=1999-12-31 06:00:00>
>>> from astroplan import Observer
>>> import datetime
>>> import pytz
>>> subaru = Observer.at_site("Subaru", timezone="US/Hawaii")
>>> hi_date_time = datetime.datetime(2005, 6, 21, 20, 0, 0, 0)
>>> subaru.datetime_to_astropy_time(hi_date_time)
<Time object: scale='utc' format='datetime' value=2005-06-22 06:00:00>
>>> utc_date_time = datetime.datetime(2005, 6, 22, 6, 0, 0, 0,
... tzinfo=pytz.timezone("UTC"))
>>> subaru.datetime_to_astropy_time(utc_date_time)
<Time object: scale='utc' format='datetime' value=2005-06-22 06:00:00>
"""
if hasattr(date_time, '__iter__'):
return Time([self.datetime_to_astropy_time(t) for t in date_time])
# For timezone-naive datetimes, assign local timezone
if date_time.tzinfo is None:
date_time = self.timezone.localize(date_time)
return Time(date_time, location=self.location)
def _is_broadcastable(self, shp1, shp2):
"""Test if two shape tuples are broadcastable"""
if shp1 == shp2:
return True
for a, b in zip(shp1[::-1], shp2[::-1]):
if a == 1 or b == 1 or a == b:
pass
else:
return False
return True
def _preprocess_inputs(self, time, target=None, grid_times_targets=False):
"""
Preprocess time and target inputs
This routine takes the inputs for time and target and attempts to
return a single `~astropy.time.Time` and `~astropy.coordinates.SkyCoord`
for each argument, which may be non-scalar if necessary.
time : `~astropy.time.Time` or other (see below)
The time(s) to use in the calculation. It can be anything that
`~astropy.time.Time` will accept (including a `~astropy.time.Time`
object)
target : `~astroplan.FixedTarget`, `~astropy.coordinates.SkyCoord`, or list
The target(s) to use in the calculation.
grid_times_targets: bool
If True, the target object will have extra dimensions packed onto
the end, so that calculations with M targets and N times will
return an (M, N) shaped result. Otherwise, we rely on broadcasting
the shapes together using standard numpy rules. Useful for grid
searches for rise/set times etc.
"""
# make sure we have a non-scalar time
if not isinstance(time, Time):
time = Time(time)
if target is None:
return time, None
# convert any kind of target argument to non-scalar SkyCoord
target = get_skycoord(target)
if grid_times_targets:
if target.isscalar:
# ensure we have a (1, 1) shape coord
target = SkyCoord(np.tile(target, 1))[:, np.newaxis]
else:
while target.ndim <= time.ndim:
target = target[:, np.newaxis]
elif not self._is_broadcastable(target.shape, time.shape):
raise ValueError('Time and Target arguments cannot be broadcast '
'against each other with shapes {} and {}'
.format(time.shape, target.shape))
return time, target
[docs] def altaz(self, time, target=None, obswl=None, grid_times_targets=False):
"""
Get an `~astropy.coordinates.AltAz` frame or coordinate.
If ``target`` is None, generates an altitude/azimuth frame. Otherwise,
calculates the transformation to that frame for the requested ``target``.
Parameters
----------
time : `~astropy.time.Time` or other (see below)
The time at which the observation is taking place. Will be used as
the ``obstime`` attribute in the resulting frame or coordinate. 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)
target : `~astroplan.FixedTarget`, `~astropy.coordinates.SkyCoord`, or list (optional)
Celestial object(s) of interest. If ``target`` is `None`, returns
the `~astropy.coordinates.AltAz` frame without coordinates.
obswl : `~astropy.units.Quantity` (optional)
Wavelength of the observation used in the calculation.
grid_times_targets: bool (optional)
If True, the target object will have extra dimensions packed
onto the end, so that calculations with M targets and N times
will return an (M, N) shaped result. Otherwise, we rely on
broadcasting the shapes together using standard numpy
rules. Useful for grid searches for rise/set times etc.
Returns
-------
`~astropy.coordinates.AltAz`
If ``target`` is `None`, returns `~astropy.coordinates.AltAz`
frame. If ``target`` is not `None`, returns the ``target``
transformed to the `~astropy.coordinates.AltAz` frame.
Examples
--------
Create an instance of the `~astropy.coordinates.AltAz` frame for an
observer at Apache Point Observatory at a particular time:
>>> from astroplan import Observer
>>> from astropy.time import Time
>>> from astropy.coordinates import SkyCoord
>>> apo = Observer.at_site("APO")
>>> time = Time('2001-02-03 04:05:06')
>>> target = SkyCoord(0*u.deg, 0*u.deg)
>>> altaz_frame = apo.altaz(time)
Now transform the target's coordinates to the alt/az frame:
>>> target_altaz = target.transform_to(altaz_frame) # doctest: +SKIP
Alternatively, construct an alt/az frame and transform the target to
that frame all in one step:
>>> target_altaz = apo.altaz(time, target) # doctest: +SKIP
"""
if target is not None:
time, target = self._preprocess_inputs(time, target, grid_times_targets)
altaz_frame = AltAz(location=self.location, obstime=time,
pressure=self.pressure, obswl=obswl,
temperature=self.temperature,
relative_humidity=self.relative_humidity)
if target is None:
# Return just the frame
return altaz_frame
else:
return target.transform_to(altaz_frame)
[docs] def parallactic_angle(self, time, target, grid_times_targets=False):
"""
Calculate the parallactic angle.
Parameters
----------
time : `~astropy.time.Time`
Observation time.
target : `~astroplan.FixedTarget` or `~astropy.coordinates.SkyCoord` or list
Target celestial object(s).
grid_times_targets: bool
If True, the target object will have extra dimensions packed onto
the end, so that calculations with M targets and N times will
return an (M, N) shaped result. Otherwise, we rely on broadcasting
the shapes together using standard numpy rules.
Returns
-------
`~astropy.coordinates.Angle`
Parallactic angle.
Notes
-----
The parallactic angle is the angle between the great circle that
intersects a celestial object and the zenith, and the object's hour
circle [1]_.
.. [1] https://en.wikipedia.org/wiki/Parallactic_angle
"""
time, coordinate = self._preprocess_inputs(time, target, grid_times_targets)
# Eqn (14.1) of Meeus' Astronomical Algorithms
LST = time.sidereal_time('mean', longitude=self.location.lon)
H = (LST - coordinate.ra).radian
q = np.arctan2(np.sin(H),
(np.tan(self.location.lat.radian) *
np.cos(coordinate.dec.radian) -
np.sin(coordinate.dec.radian)*np.cos(H)))*u.rad
return Angle(q)
# Sun-related methods.
@u.quantity_input(horizon=u.deg)
def _horiz_cross(self, t, alt, rise_set, horizon=0*u.degree):
"""
Find time ``t`` when values in array ``a`` go from
negative to positive or positive to negative (exclude endpoints)
``return_limits`` will return nearest times to zero-crossing.
Parameters
----------
t : `~astropy.time.Time`
Grid of N times, any shape. Search grid along first axis, e.g (N, ...)
alt : `~astropy.units.Quantity`
Grid of altitudes
Depending on broadcasting we either have ndim >=3 and
M targets along first axis, e.g (M, N, ...), or
ndim = 2 and targets/times in last axis
rise_set : {"rising", "setting"}
Calculate either rising or setting across the horizon
horizon : float
Number of degrees above/below actual horizon to use
for calculating rise/set times (i.e.,
-6 deg horizon = civil twilight, etc.)
Returns
-------
Returns the lower and upper limits on the time and altitudes
of the horizon crossing. The altitude limits have shape (M, ...) and the
time limits have shape (...). These arrays aresuitable for interpolation
to find the horizon crossing time.
"""
# handle different cases by enforcing standard shapes on
# the altitude grid
finesse_time_indexes = False
if alt.ndim == 1:
raise ValueError('Must supply more at least a 2D grid of altitudes')
elif alt.ndim == 2:
# TODO: this test for ndim=2 doesn't work. if times is e.g (2,5)
# then alt will have ndim=3, but shape (100, 2, 5) so grid
# is in first index...
ntargets = alt.shape[1]
ngrid = alt.shape[0]
unit = alt.unit
alt = np.broadcast_to(alt, (ntargets, ngrid, ntargets)).T
alt = alt*unit
extra_dimension_added = True
if t.shape[1] == 1:
finesse_time_indexes = True
else:
extra_dimension_added = False
output_shape = (alt.shape[0],) + alt.shape[2:]
if rise_set == 'rising':
# Find index where altitude goes from below to above horizon
condition = ((alt[:, :-1, ...] < horizon) *
(alt[:, 1:, ...] > horizon))
elif rise_set == 'setting':
# Find index where altitude goes from above to below horizon
condition = ((alt[:, :-1, ...] > horizon) *
(alt[:, 1:, ...] < horizon))
noncrossing_indices = np.sum(condition, axis=1, dtype=np.intp) < 1
alt_lims1 = u.Quantity(np.zeros(output_shape), unit=u.deg)
alt_lims2 = u.Quantity(np.zeros(output_shape), unit=u.deg)
jd_lims1 = np.zeros(output_shape)
jd_lims2 = np.zeros(output_shape)
if np.any(noncrossing_indices):
for target_index in set(np.where(noncrossing_indices)[0]):
warnmsg = ('Target with index {} does not cross horizon={} '
'within 24 hours'.format(target_index, horizon))
if (alt[target_index, ...] > horizon).all():
warnings.warn(warnmsg, TargetAlwaysUpWarning)
else:
warnings.warn(warnmsg, TargetNeverUpWarning)
alt_lims1[np.nonzero(noncrossing_indices)] = np.nan
alt_lims2[np.nonzero(noncrossing_indices)] = np.nan
jd_lims1[np.nonzero(noncrossing_indices)] = np.nan
jd_lims2[np.nonzero(noncrossing_indices)] = np.nan
before_indices = np.array(np.nonzero(condition))
# we want to add an vector like (0, 1, ...) to get after indices
after_indices = before_indices.copy()
after_indices[1, :] += 1
al1 = alt[tuple(before_indices)]
al2 = alt[tuple(after_indices)]
# slice the time in the same way, but delete the object index
before_time_index_tuple = np.delete(before_indices, 0, 0)
after_time_index_tuple = np.delete(after_indices, 0, 0)
if finesse_time_indexes:
before_time_index_tuple[1:] = 0
after_time_index_tuple[1:] = 0
tl1 = t[tuple(before_time_index_tuple)]
tl2 = t[tuple(after_time_index_tuple)]
alt_lims1[tuple(np.delete(before_indices, 1, 0))] = al1
alt_lims2[tuple(np.delete(before_indices, 1, 0))] = al2
jd_lims1[tuple(np.delete(before_indices, 1, 0))] = tl1.utc.jd
jd_lims2[tuple(np.delete(before_indices, 1, 0))] = tl2.utc.jd
if extra_dimension_added:
return (alt_lims1.diagonal(), alt_lims2.diagonal(),
jd_lims1.diagonal(), jd_lims2.diagonal())
else:
return alt_lims1, alt_lims2, jd_lims1, jd_lims2
@u.quantity_input(horizon=u.deg)
def _two_point_interp(self, jd_before, jd_after,
alt_before, alt_after, horizon=0*u.deg):
"""
Do linear interpolation between two ``altitudes`` at
two ``times`` to determine the time where the altitude
goes through zero.
Parameters
----------
jd_before : `float`
JD(UTC) before crossing event
jd_after : `float`
JD(UTC) after crossing event
alt_before : `~astropy.units.Quantity`
altitude before crossing event
alt_after : `~astropy.units.Quantity`
altitude after crossing event
horizon : `~astropy.units.Quantity`
Solve for the time when the altitude is equal to
reference_alt.
Returns
-------
t : `~astropy.time.Time`
Time when target crosses the horizon
"""
# Approximate the horizon-crossing time:
slope = (alt_after-alt_before) / ((jd_after - jd_before) * u.d)
crossing_jd = (jd_after * u.d - ((alt_after - horizon) / slope))
# TODO: edit after https://github.com/astropy/astropy/issues/9612 has
# been addressed.
# Determine whether or not there are NaNs in the crossing_jd array which
# represent computations where no horizon crossing was found:
nans = np.isnan(crossing_jd)
# If there are, set them equal to zero, rather than np.nan
crossing_jd[nans] = 0
times = Time(crossing_jd, format='jd')
# Create a Time object with masked out times where there were NaNs
times[nans] = np.ma.masked
return np.squeeze(times)
def _altitude_trig(self, LST, target, grid_times_targets=False):
"""
Calculate the altitude of ``target`` at local sidereal times ``LST``.
This method provides a factor of ~3 speed up over calling `altaz`, and
inherently does *not* take the atmosphere into account.
Parameters
----------
LST : `~astropy.time.Time`
Local sidereal times (array)
target : {`~astropy.coordinates.SkyCoord`, `FixedTarget`} or similar
Target celestial object's coordinates.
grid_times_targets: bool
If True, the target object will have extra dimensions packed onto
the end, so that calculations with M targets and N times will
return an (M, N) shaped result. Otherwise, we rely on broadcasting
the shapes together using standard numpy rules.
Returns
-------
alt : `~astropy.unit.Quantity`
Array of altitudes
"""
LST, target = self._preprocess_inputs(LST, target, grid_times_targets)
alt = np.arcsin(np.sin(self.location.lat.radian) *
np.sin(target.dec) +
np.cos(self.location.lat.radian) *
np.cos(target.dec) *
np.cos(LST.radian - target.ra.radian))
return alt
def _calc_riseset(self, time, target, prev_next, rise_set, horizon,
n_grid_points=150, grid_times_targets=False):
"""
Time at next rise/set of ``target``.
Parameters
----------
time : `~astropy.time.Time` or other (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)
target : `~astropy.coordinates.SkyCoord`
Position of target or multiple positions of that target
at multiple times (if target moves, like the Sun)
prev_next : str - either 'previous' or 'next'
Test next rise/set or previous rise/set
rise_set : str - either 'rising' or 'setting'
Compute prev/next rise or prev/next set
horizon : `~astropy.units.Quantity`
Degrees above/below actual horizon to use
for calculating rise/set times (i.e.,
-6 deg horizon = civil twilight, etc.)
n_grid_points : int (optional)
Number of altitudes to compute when searching for
rise or set.
grid_times_targets: bool
If True, the target object will have extra dimensions packed onto
the end, so that calculations with M targets and N times will
return an (M, N) shaped result. Otherwise, we rely on broadcasting
the shapes together using standard numpy rules.
Returns
-------
ret1 : `~astropy.time.Time`
Time of rise/set
"""
if not isinstance(time, Time):
time = Time(time)
if prev_next == 'next':
start = 0
end = (1 + (target.approx_sidereal_drift.to(u.day).value
if hasattr(target, 'approx_sidereal_drift') else 0))
else:
start = (-1 - (target.approx_sidereal_drift.to(u.day).value
if hasattr(target, 'approx_sidereal_drift') else 0))
end = 0
times = _generate_24hr_grid(time, start, end, n_grid_points)
if target is MoonFlag:
altaz = self.altaz(times, get_moon(times, location=self.location),
grid_times_targets=grid_times_targets)
elif target is SunFlag:
altaz = self.altaz(times, get_sun(times),
grid_times_targets=grid_times_targets)
else:
altaz = self.altaz(times, target,
grid_times_targets=grid_times_targets)
altitudes = altaz.alt
al1, al2, jd1, jd2 = self._horiz_cross(times, altitudes, rise_set,
horizon)
return self._two_point_interp(jd1, jd2, al1, al2,
horizon=horizon)
def _calc_transit(self, time, target, prev_next, antitransit=False,
n_grid_points=150, grid_times_targets=False):
"""
Time at next transit of the meridian of `target`.
Parameters
----------
time : `~astropy.time.Time` or other (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)
target : `~astropy.coordinates.SkyCoord`
Position of target or multiple positions of that target
at multiple times (if target moves, like the Sun)
prev_next : str - either 'previous' or 'next'
Test next rise/set or previous rise/set
antitransit : bool
Toggle compute antitransit (below horizon, equivalent to midnight
for the Sun)
n_grid_points : int (optional)
Number of altitudes to compute when searching for
rise or set.
grid_times_targets: bool
If True, the target object will have extra dimensions packed onto
the end, so that calculations with M targets and N times will
return an (M, N) shaped result. Otherwise, we rely on broadcasting
the shapes together using standard numpy rules.
Returns
-------
ret1 : `~astropy.time.Time`
Time of transit/antitransit
"""
# TODO FIX BROADCASTING HERE
if not isinstance(time, Time):
time = Time(time)
if prev_next == 'next':
times = _generate_24hr_grid(time, 0, 1, n_grid_points,
for_deriv=True)
else:
times = _generate_24hr_grid(time, -1, 0, n_grid_points,
for_deriv=True)
# The derivative of the altitude with respect to time is increasing
# from negative to positive values at the anti-transit of the meridian
if antitransit:
rise_set = 'rising'
else:
rise_set = 'setting'
altaz = self.altaz(times, target, grid_times_targets=grid_times_targets)
altitudes = altaz.alt
if altitudes.ndim > 2:
# shape is (M, N, ...) where M is targets and N is grid
d_altitudes = altitudes.diff(axis=1)
else:
# shape is (N, M) where M is targets and N is grid
d_altitudes = altitudes.diff(axis=0)
dt = Time((times.jd[1:] + times.jd[:-1])/2, format='jd')
horizon = 0*u.degree # Find when derivative passes through zero
al1, al2, jd1, jd2 = self._horiz_cross(dt, d_altitudes,
rise_set, horizon)
return self._two_point_interp(jd1, jd2, al1, al2,
horizon=horizon)
def _determine_which_event(self, function, args_dict):
"""
Run through the next/previous/nearest permutations of the solutions
to `function(time, ...)`, and return the previous/next/nearest one
specified by the args stored in args_dict.
"""
time = args_dict.pop('time', None)
target = args_dict.pop('target', None)
which = args_dict.pop('which', None)
horizon = args_dict.pop('horizon', None)
rise_set = args_dict.pop('rise_set', None)
antitransit = args_dict.pop('antitransit', None)
grid_times_targets = args_dict.pop('grid_times_targets', False)
n_grid_points = args_dict.pop('n_grid_points', 150)
# Assemble arguments for function, depending on the function.
if function == self._calc_riseset:
def event_function(w):
return function(time, target, w, rise_set, horizon,
grid_times_targets=grid_times_targets,
n_grid_points=n_grid_points)
elif function == self._calc_transit:
def event_function(w):
return function(time, target, w, antitransit=antitransit,
grid_times_targets=grid_times_targets,
n_grid_points=n_grid_points)
else:
raise ValueError('Function {} not supported in '
'_determine_which_event.'.format(function))
if not isinstance(time, Time):
time = Time(time)
if which == 'next' or which == 'nearest':
next_event = event_function('next')
if which == 'next':
return next_event
if which == 'previous' or which == 'nearest':
previous_event = event_function('previous')
if which == 'previous':
return previous_event
if which == 'nearest':
# Use some hacks to handle the non-rising/non-setting cases
try:
mask = abs(time - previous_event) < abs(time - next_event)
except TypeError:
# encountered if time is scalar & nan
return next_event
ma = np.where(mask, previous_event.utc.jd, next_event.utc.jd)
# HACK: Time objects cannot be initiated w/NaN, so we first
# make them zero, then change them to NaN
not_finite = ~np.isfinite(ma)
ma[not_finite] = 0
tm = Time(ma, format='jd')
tm[not_finite] = np.nan
return tm
raise ValueError('"which" kwarg must be "next", "previous" or '
'"nearest".')
[docs] @u.quantity_input(horizon=u.deg)
def target_rise_time(self, time, target, which='nearest',
horizon=0*u.degree, grid_times_targets=False, n_grid_points=150):
"""
Calculate rise time.
Compute time of the next/previous/nearest rise of the ``target``
object, where "rise" is defined as the time when the ``target``
transitions from altitudes below the ``horizon`` to above the
``horizon``.
Parameters
----------
time : `~astropy.time.Time` or other (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)
target : `~astropy.coordinates.SkyCoord`, `~astroplan.FixedTarget`, or list
Target celestial object(s)
which : {'next', 'previous', 'nearest'}
Choose which sunrise relative to the present ``time`` would you
like to calculate
horizon : `~astropy.units.Quantity` (optional), default = zero degrees
Degrees above/below actual horizon to use
for calculating rise/set times (i.e.,
-6 deg horizon = civil twilight, etc.)
grid_times_targets: bool
If True, the target object will have extra dimensions packed
onto the end, so that calculations with M targets and N times
will return an (M, N) shaped result. Otherwise, we rely on
broadcasting the shapes together using standard numpy rules.
n_grid_points : int (optional)
The number of grid points on which to search for the horizon
crossings of the target over a 24 hour period, default is 150 which
yields rise time precisions better than one minute.
Returns
-------
`~astropy.time.Time`
Rise time of target
Examples
--------
Calculate the rise time of Rigel at Keck Observatory:
>>> from astroplan import Observer, FixedTarget
>>> from astropy.time import Time
>>> time = Time("2001-02-03 04:05:06")
>>> target = FixedTarget.from_name("Rigel")
>>> keck = Observer.at_site("Keck")
>>> rigel_rise_time = keck.target_rise_time(time, target, which="next") # doctest: +SKIP
>>> print("ISO: {0.iso}, JD: {0.jd}".format(rigel_rise_time)) # doctest: +SKIP
ISO: 2001-02-04 00:51:23.330, JD: 2451944.53569
"""
return self._determine_which_event(self._calc_riseset,
dict(time=time, target=target,
which=which, rise_set='rising',
horizon=horizon,
n_grid_points=n_grid_points,
grid_times_targets=grid_times_targets))
[docs] @u.quantity_input(horizon=u.deg)
def target_set_time(self, time, target, which='nearest', horizon=0*u.degree,
grid_times_targets=False, n_grid_points=150):
"""
Calculate set time.
Compute time of the next/previous/nearest set of ``target``, where
"set" is defined as when the ``target`` transitions from altitudes
above ``horizon`` to below ``horizon``.
Parameters
----------
time : `~astropy.time.Time` or other (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)
target : `~astropy.coordinates.SkyCoord`, `~astroplan.FixedTarget`, or list
Target celestial object(s)
which : {'next', 'previous', 'nearest'}
Choose which sunset relative to the present ``time`` would you
like to calculate
horizon : `~astropy.units.Quantity` (optional), default = zero degrees
Degrees above/below actual horizon to use
for calculating rise/set times (i.e.,
-6 deg horizon = civil twilight, etc.)
grid_times_targets: bool
If True, the target object will have extra dimensions packed
onto the end, so that calculations with M targets and N times
will return an (M, N) shaped result. Otherwise, we rely on
broadcasting the shapes together using standard numpy rules.
n_grid_points : int (optional)
The number of grid points on which to search for the horizon
crossings of the target over a 24 hour period, default is 150 which
yields set time precisions better than one minute.
Returns
-------
`~astropy.time.Time`
Set time of target.
Examples
--------
Calculate the set time of Rigel at Keck Observatory:
>>> from astroplan import Observer, FixedTarget
>>> from astropy.time import Time
>>> time = Time("2001-02-03 04:05:06")
>>> target = FixedTarget.from_name("Rigel")
>>> keck = Observer.at_site("Keck")
>>> rigel_set_time = keck.target_set_time(time, target, which="next") # doctest: +SKIP
>>> print("ISO: {0.iso}, JD: {0.jd}".format(rigel_set_time)) # doctest: +SKIP
ISO: 2001-02-03 12:29:34.768, JD: 2451944.02054
"""
return self._determine_which_event(self._calc_riseset,
dict(time=time, target=target,
which=which,
rise_set='setting',
horizon=horizon,
n_grid_points=n_grid_points,
grid_times_targets=grid_times_targets))
[docs] def target_meridian_transit_time(self, time, target, which='nearest',
grid_times_targets=False, n_grid_points=150):
"""
Calculate time at the transit of the meridian.
Compute time of the next/previous/nearest transit of the ``target``
object.
Parameters
----------
time : `~astropy.time.Time` or other (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)
target : `~astropy.coordinates.SkyCoord`, `~astroplan.FixedTarget`, or list
Target celestial object(s)
which : {'next', 'previous', 'nearest'}
Choose which sunrise relative to the present ``time`` would you
like to calculate
grid_times_targets: bool
If True, the target object will have extra dimensions packed
onto the end, so that calculations with M targets and N times
will return an (M, N) shaped result. Otherwise, we rely on
broadcasting the shapes together using standard numpy rules.
n_grid_points : int (optional)
The number of grid points on which to search for the horizon
crossings of the target over a 24 hour period, default is 150 which
yields rise time precisions better than one minute.
Returns
-------
`~astropy.time.Time`
Transit time of target
Examples
--------
Calculate the meridian transit time of Rigel at Keck Observatory:
>>> from astroplan import Observer, FixedTarget
>>> from astropy.time import Time
>>> time = Time("2001-02-03 04:05:06")
>>> target = FixedTarget.from_name("Rigel")
>>> keck = Observer.at_site("Keck")
>>> rigel_transit_time = keck.target_meridian_transit_time(time, target,
... which="next") # doctest: +SKIP
>>> print("ISO: {0.iso}, JD: {0.jd}".format(rigel_transit_time)) # doctest: +SKIP
ISO: 2001-02-03 06:42:26.863, JD: 2451943.77948
"""
return self._determine_which_event(self._calc_transit,
dict(time=time, target=target,
which=which,
n_grid_points=n_grid_points,
rise_set='setting',
grid_times_targets=grid_times_targets))
[docs] def target_meridian_antitransit_time(self, time, target, which='nearest',
grid_times_targets=False, n_grid_points=150):
"""
Calculate time at the antitransit of the meridian.
Compute time of the next/previous/nearest antitransit of the ``target``
object.
Parameters
----------
time : `~astropy.time.Time` or other (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)
target : `~astropy.coordinates.SkyCoord`, `~astroplan.FixedTarget`, or list
Target celestial object(s)
which : {'next', 'previous', 'nearest'}
Choose which sunrise relative to the present ``time`` would you
like to calculate
grid_times_targets : bool
If True, the target object will have extra dimensions packed onto the end,
so that calculations with M targets and N times will return an (M, N)
shaped result. Otherwise, we rely on broadcasting the shapes together
using standard numpy rules.
n_grid_points : int (optional)
The number of grid points on which to search for the horizon
crossings of the target over a 24 hour period, default is 150 which
yields rise time precisions better than one minute.
Returns
-------
`~astropy.time.Time`
Antitransit time of target
Examples
--------
Calculate the meridian anti-transit time of Rigel at Keck Observatory:
>>> from astroplan import Observer, FixedTarget
>>> from astropy.time import Time
>>> time = Time("2001-02-03 04:05:06")
>>> target = FixedTarget.from_name("Rigel")
>>> keck = Observer.at_site("Keck")
>>> rigel_antitransit_time = keck.target_meridian_antitransit_time(
... time, target, which="next") # doctest: +SKIP
>>> print("ISO: {0.iso}, JD: {0.jd}".format(rigel_antitransit_time)) # doctest: +SKIP
ISO: 2001-02-03 18:40:29.761, JD: 2451944.27812
"""
return self._determine_which_event(self._calc_transit,
dict(time=time, target=target,
which=which, antitransit=True,
rise_set='setting',
n_grid_points=n_grid_points,
grid_times_targets=grid_times_targets))
[docs] @u.quantity_input(horizon=u.deg)
def sun_rise_time(self, time, which='nearest', horizon=0*u.degree, n_grid_points=150):
"""
Time of sunrise.
Compute time of the next/previous/nearest sunrise, where
sunrise is defined as when the Sun transitions from altitudes
below ``horizon`` to above ``horizon``.
Parameters
----------
time : `~astropy.time.Time` or other (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)
which : {'next', 'previous', 'nearest'}
Choose which sunrise relative to the present ``time`` would you
like to calculate.
horizon : `~astropy.units.Quantity` (optional), default = zero degrees
Degrees above/below actual horizon to use
for calculating rise/set times (i.e.,
-6 deg horizon = civil twilight, etc.)
n_grid_points : int (optional)
The number of grid points on which to search for the horizon
crossings of the target over a 24 hour period, default is 150 which
yields rise time precisions better than one minute.
Returns
-------
`~astropy.time.Time`
Time of sunrise
Examples
--------
Calculate the time of the previous sunrise at Apache Point Observatory:
>>> from astroplan import Observer
>>> from astropy.time import Time
>>> apo = Observer.at_site("APO")
>>> time = Time('2001-02-03 04:05:06')
>>> sun_rise = apo.sun_rise_time(time, which="previous") # doctest: +SKIP
>>> print("ISO: {0.iso}, JD: {0.jd}".format(sun_rise)) # doctest: +SKIP
ISO: 2001-02-02 14:02:50.554, JD: 2451943.08531
"""
return self.target_rise_time(time, get_sun(time), which, horizon,
n_grid_points=n_grid_points)
[docs] @u.quantity_input(horizon=u.deg)
def sun_set_time(self, time, which='nearest', horizon=0*u.degree, n_grid_points=150):
"""
Time of sunset.
Compute time of the next/previous/nearest sunset, where
sunset is defined as when the Sun transitions from altitudes
below ``horizon`` to above ``horizon``.
Parameters
----------
time : `~astropy.time.Time` or other (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)
which : {'next', 'previous', 'nearest'}
Choose which sunset relative to the present ``time`` would you
like to calculate
horizon : `~astropy.units.Quantity` (optional), default = zero degrees
Degrees above/below actual horizon to use
for calculating rise/set times (i.e.,
-6 deg horizon = civil twilight, etc.)
n_grid_points : int (optional)
The number of grid points on which to search for the horizon
crossings of the target over a 24 hour period, default is 150 which
yields set time precisions better than one minute.
Returns
-------
`~astropy.time.Time`
Time of sunset
Examples
--------
Calculate the time of the next sunset at Apache Point Observatory:
>>> from astroplan import Observer
>>> from astropy.time import Time
>>> apo = Observer.at_site("APO")
>>> time = Time('2001-02-03 04:05:06')
>>> sun_set = apo.sun_set_time(time, which="next") # doctest: +SKIP
>>> print("ISO: {0.iso}, JD: {0.jd}".format(sun_set)) # doctest: +SKIP
ISO: 2001-02-04 00:35:42.102, JD: 2451944.52479
"""
return self.target_set_time(time, get_sun(time), which, horizon,
n_grid_points=n_grid_points)
[docs] def noon(self, time, which='nearest', n_grid_points=150):
"""
Time at solar noon.
Parameters
----------
time : `~astropy.time.Time` or other (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)
which : {'next', 'previous', 'nearest'}
Choose which noon relative to the present ``time`` would you
like to calculate
n_grid_points : int (optional)
The number of grid points on which to search for the horizon
crossings of the target over a 24 hour period, default is 150 which
yields noon time precisions better than one minute.
Returns
-------
`~astropy.time.Time`
Time at solar noon
"""
return self.target_meridian_transit_time(time, get_sun(time), which,
n_grid_points=n_grid_points)
[docs] def midnight(self, time, which='nearest', n_grid_points=150):
"""
Time at solar midnight.
Parameters
----------
time : `~astropy.time.Time` or other (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)
which : {'next', 'previous', 'nearest'}
Choose which noon relative to the present ``time`` would you
like to calculate
n_grid_points : int (optional)
The number of grid points on which to search for the horizon
crossings of the target over a 24 hour period, default is 150 which
yields midnight time precisions better than one minute.
Returns
-------
`~astropy.time.Time`
Time at solar midnight
"""
return self.target_meridian_antitransit_time(time, get_sun(time), which,
n_grid_points=n_grid_points)
# Twilight convenience functions
[docs] def twilight_evening_astronomical(self, time, which='nearest', n_grid_points=150):
"""
Time at evening astronomical (-18 degree) twilight.
Parameters
----------
time : `~astropy.time.Time` or other (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)
which : {'next', 'previous', 'nearest'}
Choose which twilight relative to the present ``time`` would you
like to calculate. Default is nearest.
n_grid_points : int (optional)
The number of grid points on which to search for the horizon
crossings of the target over a 24 hour period, default is 150 which
yields twilight time precisions better than one minute.
Returns
-------
`~astropy.time.Time`
Time of twilight
"""
return self.sun_set_time(time, which, horizon=-18*u.degree,
n_grid_points=n_grid_points)
[docs] def twilight_evening_nautical(self, time, which='nearest', n_grid_points=150):
"""
Time at evening nautical (-12 degree) twilight.
Parameters
----------
time : `~astropy.time.Time` or other (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)
which : {'next', 'previous', 'nearest'}
Choose which twilight relative to the present ``time`` would you
like to calculate. Default is nearest.
n_grid_points : int (optional)
The number of grid points on which to search for the horizon
crossings of the target over a 24 hour period, default is 150 which
yields twilight time precisions better than one minute.
Returns
-------
`~astropy.time.Time`
Time of twilight
"""
return self.sun_set_time(time, which, horizon=-12*u.degree,
n_grid_points=n_grid_points)
[docs] def twilight_evening_civil(self, time, which='nearest', n_grid_points=150):
"""
Time at evening civil (-6 degree) twilight.
Parameters
----------
time : `~astropy.time.Time` or other (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)
which : {'next', 'previous', 'nearest'}
Choose which twilight relative to the present ``time`` would you
like to calculate. Default is nearest.
n_grid_points : int (optional)
The number of grid points on which to search for the horizon
crossings of the target over a 24 hour period, default is 150 which
yields twilight time precisions better than one minute.
Returns
-------
`~astropy.time.Time`
Time of twilight
"""
return self.sun_set_time(time, which, horizon=-6*u.degree,
n_grid_points=n_grid_points)
[docs] def twilight_morning_astronomical(self, time, which='nearest', n_grid_points=150):
"""
Time at morning astronomical (-18 degree) twilight.
Parameters
----------
time : `~astropy.time.Time` or other (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)
which : {'next', 'previous', 'nearest'}
Choose which twilight relative to the present ``time`` would you
like to calculate
n_grid_points : int (optional)
The number of grid points on which to search for the horizon
crossings of the target over a 24 hour period, default is 150 which
yields twilight time precisions better than one minute.
Returns
-------
`~astropy.time.Time`
Time of twilight
"""
return self.sun_rise_time(time, which, horizon=-18*u.degree,
n_grid_points=n_grid_points)
[docs] def twilight_morning_nautical(self, time, which='nearest', n_grid_points=150):
"""
Time at morning nautical (-12 degree) twilight.
Parameters
----------
time : `~astropy.time.Time` or other (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)
which : {'next', 'previous', 'nearest'}
Choose which twilight relative to the present ``time`` would you
like to calculate. Default is nearest.
n_grid_points : int (optional)
The number of grid points on which to search for the horizon
crossings of the target over a 24 hour period, default is 150 which
yields twilight time precisions better than one minute.
Returns
-------
`~astropy.time.Time`
Time of twilight
"""
return self.sun_rise_time(time, which, horizon=-12*u.degree,
n_grid_points=n_grid_points)
[docs] def twilight_morning_civil(self, time, which='nearest', n_grid_points=150):
"""
Time at morning civil (-6 degree) twilight.
Parameters
----------
time : `~astropy.time.Time` or other (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)
which : {'next', 'previous', 'nearest'}
Choose which twilight relative to the present ``time`` would you
like to calculate. Default is nearest.
n_grid_points : int (optional)
The number of grid points on which to search for the horizon
crossings of the target over a 24 hour period, default is 150 which
yields twilight time precisions better than one minute.
Returns
-------
`~astropy.time.Time`
Time of sunset
"""
return self.sun_rise_time(time, which, horizon=-6*u.degree,
n_grid_points=n_grid_points)
# Moon-related methods.
[docs] def moon_rise_time(self, time, which='nearest', horizon=0*u.deg, n_grid_points=150):
"""
Returns the local moon rise time.
Compute time of the next/previous/nearest moon rise, where
moon rise is defined as the time when the moon transitions from
altitudes below ``horizon`` to above ``horizon``.
Parameters
----------
time : `~astropy.time.Time` or other (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)
which : {'next', 'previous', 'nearest'}
Choose which moon rise relative to the present ``time`` would you
like to calculate.
horizon : `~astropy.units.Quantity` (optional), default = zero degrees
Degrees above/below actual horizon to use
for calculating rise/set times (i.e.,
-6 deg horizon = civil twilight, etc.)
n_grid_points : int (optional)
The number of grid points on which to search for the horizon
crossings of the target over a 24 hour period, default is 150 which
yields rise time precisions better than one minute.
"""
return self.target_rise_time(time, MoonFlag, which, horizon,
n_grid_points=n_grid_points)
[docs] def moon_set_time(self, time, which='nearest', horizon=0*u.deg, n_grid_points=150):
"""
Returns the local moon set time.
Compute time of the next/previous/nearest moon set, where
moon set is defined as the time when the moon transitions from
altitudes below ``horizon`` to above ``horizon``.
Parameters
----------
time : `~astropy.time.Time` or other (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)
which : {'next', 'previous', 'nearest'}
Choose which moon set relative to the present ``time`` would you
like to calculate.
horizon : `~astropy.units.Quantity` (optional), default = zero degrees
Degrees above/below actual horizon to use
for calculating set/set times (i.e.,
-6 deg horizon = civil twilight, etc.)
n_grid_points : int (optional)
The number of grid points on which to search for the horizon
crossings of the target over a 24 hour period, default is 150 which
yields set time precisions better than one minute.
"""
return self.target_set_time(time, MoonFlag, which, horizon,
n_grid_points=n_grid_points)
[docs] def moon_illumination(self, time):
"""
Calculate the illuminated fraction of the moon.
Parameters
----------
time : `~astropy.time.Time` or other (see below)
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).
Returns
-------
float
Fraction of lunar surface illuminated
Examples
--------
How much of the lunar surface is illuminated at 2015-08-29 18:35 UTC,
which we happen to know is the time of a full moon?
>>> from astroplan import Observer
>>> from astropy.time import Time
>>> apo = Observer.at_site("APO")
>>> time = Time("2015-08-29 18:35")
>>> apo.moon_illumination(time) # doctest: +SKIP
array([ 0.99972487])
"""
if not isinstance(time, Time):
time = Time(time)
return moon_illumination(time)
[docs] def moon_phase(self, time=None):
"""
Calculate lunar orbital phase.
For example, phase=pi is "new", phase=0 is "full".
Parameters
----------
time : `~astropy.time.Time` or other (see below)
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).
Returns
-------
moon_phase_angle : float
Orbital phase angle of the moon where pi corresponds to new moon,
zero corresponds to full moon.
Examples
--------
Calculate the phase of the moon at 2015-08-29 18:35 UTC. Near zero
radians corresponds to a nearly full moon.
>>> from astroplan import Observer
>>> from astropy.time import Time
>>> apo = Observer.at_site('APO')
>>> time = Time('2015-08-29 18:35')
>>> apo.moon_phase(time) # doctest: +SKIP
<Quantity [ 0.03317537] rad>
"""
if time is not None and not isinstance(time, Time):
time = Time(time)
return moon_phase_angle(time)
[docs] def moon_altaz(self, time, ephemeris=None):
"""
Returns the position of the moon in alt/az.
Parameters
----------
time : `~astropy.time.Time` or other (see below)
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).
ephemeris : str, optional
Ephemeris to use. If not given, use the one set with
``astropy.coordinates.solar_system_ephemeris.set`` (which is
set to 'builtin' by default).
Returns
-------
altaz : `~astropy.coordinates.SkyCoord`
Position of the moon transformed to altitude and azimuth
Examples
--------
Calculate the altitude and azimuth of the moon at Apache Point
Observatory:
>>> from astroplan import Observer
>>> from astropy.time import Time
>>> apo = Observer.at_site("APO")
>>> time = Time("2015-08-29 18:35")
>>> altaz_moon = apo.moon_altaz(time) # doctest: +SKIP
>>> print("alt: {0.alt}, az: {0.az}".format(altaz_moon)) # doctest: +SKIP
alt: -63.72706397691421 deg, az: 345.3640380598265 deg
"""
if not isinstance(time, Time):
time = Time(time)
moon = get_moon(time, location=self.location, ephemeris=ephemeris)
return self.altaz(time, moon)
[docs] def sun_altaz(self, time):
"""
Returns the position of the Sun in alt/az.
Parameters
----------
time : `~astropy.time.Time` or other (see below)
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).
ephemeris : str, optional
Ephemeris to use. If not given, use the one set with
``astropy.coordinates.solar_system_ephemeris.set`` (which is
set to 'builtin' by default).
Returns
-------
altaz : `~astropy.coordinates.SkyCoord`
Position of the sun transformed to altitude and azimuth
"""
if not isinstance(time, Time):
time = Time(time)
sun = get_sun(time)
return self.altaz(time, sun)
[docs] @u.quantity_input(horizon=u.deg)
def target_is_up(self, time, target, horizon=0*u.degree,
return_altaz=False, grid_times_targets=False):
"""
Is ``target`` above ``horizon`` at this ``time``?
Parameters
----------
time : `~astropy.time.Time` or other (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)
target : `~astropy.coordinates.SkyCoord`, `~astroplan.FixedTarget`, or list
Target celestial object(s)
horizon : `~astropy.units.Quantity` (optional), default = zero degrees
Degrees above/below actual horizon to use
for calculating rise/set times (i.e.,
-6 deg horizon = civil twilight, etc.)
return_altaz : bool (optional)
Also return the '~astropy.coordinates.AltAz' coordinate.
grid_times_targets: bool
If True, the target object will have extra dimensions packed
onto the end, so that calculations with M targets and N times
will return an (M, N) shaped result. Otherwise, we rely on
broadcasting the shapes together using standard numpy rules.
Returns
-------
observable : boolean or np.ndarray(bool)
True if ``target`` is above ``horizon`` at ``time``, else False.
Examples
--------
Are Aldebaran and Vega above the horizon at Apache Point Observatory
at 2015-08-29 18:35 UTC?
>>> from astroplan import Observer, FixedTarget
>>> from astropy.time import Time
>>> apo = Observer.at_site("APO")
>>> time = Time("2015-08-29 18:35")
>>> aldebaran = FixedTarget.from_name("Aldebaran")
>>> vega = FixedTarget.from_name("Vega")
>>> apo.target_is_up(time, aldebaran) # doctest: +SKIP
True
>>> apo.target_is_up(time, [aldebaran, vega]) # doctest: +SKIP
array([ True, False], dtype=bool)
"""
if not isinstance(time, Time):
time = Time(time)
altaz = self.altaz(time, target, grid_times_targets=grid_times_targets)
observable = altaz.alt > horizon
if altaz.isscalar:
observable = bool(observable)
if not return_altaz:
return observable
else:
return observable, altaz
[docs] @u.quantity_input(horizon=u.deg)
def is_night(self, time, horizon=0*u.deg, obswl=None):
"""
Is the Sun below ``horizon`` at ``time``?
Parameters
----------
time : `~astropy.time.Time` or other (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)
horizon : `~astropy.units.Quantity` (optional), default = zero degrees
Degrees above/below actual horizon to use
for calculating day/night (i.e.,
-6 deg horizon = civil twilight, etc.)
obswl : `~astropy.units.Quantity` (optional)
Wavelength of the observation used in the calculation
Returns
-------
sun_below_horizon : bool or np.ndarray(bool)
`True` if sun is below ``horizon`` at ``time``, else `False`.
Examples
--------
Is it "nighttime" (i.e. is the Sun below ``horizon``) at Apache Point
Observatory at 2015-08-29 18:35 UTC?
>>> from astroplan import Observer
>>> from astropy.time import Time
>>> apo = Observer.at_site("APO")
>>> time = Time("2015-08-29 18:35")
>>> apo.is_night(time) # doctest: +SKIP
False
"""
if not isinstance(time, Time):
time = Time(time)
solar_altitude = self.altaz(time, target=get_sun(time), obswl=obswl).alt
if solar_altitude.isscalar:
return bool(solar_altitude < horizon)
else:
return solar_altitude < horizon
[docs] def local_sidereal_time(self, time, kind='apparent', model=None):
"""
Convert ``time`` to local sidereal time for observer.
This is a thin wrapper around the `~astropy.time.Time.sidereal_time`
method.
Parameters
----------
time : `~astropy.time.Time` or other (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)
kind : {'mean', 'apparent'} (optional)
Passed to the ``kind`` argument of
`~astropy.time.Time.sidereal_time`
model : str or `None`; optional
The precession/nutation model to assume - see
`~astropy.time.Time.sidereal_time` for more details.
Returns
-------
`~astropy.coordinates.Longitude`
Local sidereal time.
"""
if not isinstance(time, Time):
time = Time(time)
return time.sidereal_time(kind, longitude=self.location.lon,
model=model)
[docs] def target_hour_angle(self, time, target, grid_times_targets=False):
"""
Calculate the local hour angle of ``target`` at ``time``.
Parameters
----------
time : `~astropy.time.Time` or other (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)
target : `~astropy.coordinates.SkyCoord`, `~astroplan.FixedTarget`, or list
Target celestial object(s)
grid_times_targets: bool
If True, the target object will have extra dimensions packed
onto the end, so that calculations with M targets and N times
will return an (M, N) shaped result. Otherwise, we rely on
broadcasting the shapes together using standard numpy rules.
Returns
-------
hour_angle : `~astropy.coordinates.Angle`
The hour angle(s) of the target(s) at ``time``
"""
time, target = self._preprocess_inputs(time, target, grid_times_targets)
return Longitude(self.local_sidereal_time(time) - target.ra)
[docs] @u.quantity_input(horizon=u.degree)
def tonight(self, time=None, horizon=0 * u.degree, obswl=None):
"""
Return a time range corresponding to the nearest night
This will return a range of `~astropy.time.Time` corresponding to the
beginning and ending of the night. If in the middle of a given night,
return times from `~astropy.time.Time.now` until the nearest
`~astroplan.Observer.sun_rise_time`
Parameters
----------
time : `~astropy.time.Time` (optional), default = `~astropy.time.Time.now`
The start time for tonight, which is allowed to be arbitrary.
See description above for behavior
horizon : `~astropy.units.Quantity` (optional), default = zero degrees
Degrees above/below actual horizon to use for calculating rise/set
times (e.g., -6 deg horizon = civil twilight, etc.)
obswl : `~astropy.units.Quantity` (optional)
Wavelength of the observation used in the calculation
Returns
-------
times : `~astropy.time.Time`
A tuple of times corresponding to the start and end of current night
"""
current_time = Time.now() if time is None else time
night_mask = self.is_night(current_time, horizon=horizon, obswl=obswl)
sun_set_time = self.sun_set_time(current_time, which='next',
horizon=horizon)
start_time = np.where(night_mask, current_time, sun_set_time)
# np.where gives us a list of start Times - convert to Time object
if not isinstance(start_time, Time):
start_time = Time(start_time)
end_time = self.sun_rise_time(start_time, which='next', horizon=horizon)
return start_time, end_time