In [72]:
from __future__ import (absolute_import, division, print_function,
unicode_literals)
from astropy.coordinates import (EarthLocation, Latitude, Longitude,
SkyCoord, get_sun, AltAz)
import astropy.units as u
from astropy.time import Time
import numpy as np
################################################################################
# TODO: Temporary solution to IERS tables problems
from astropy.utils.data import download_file
from astropy.utils import iers
import datetime
iers.IERS.iers_table = iers.IERS_A.open(download_file(iers.IERS_A_URL, cache=True))
################################################################################
def _horiz_cross(t, alt, rise_set, horizon=0, unit=u.degree,
return_limits=False, return_alts=False):
'''
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 times
alt : `~astropy.units.Quantity`
Grid of altitudes
rise_set : str, either "rising" or "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.)
unit : astropy unit
Unit for `horizon`
return_limits : bool
If `True`, return the lower and upper limits on the
horizon crossing time. If `False`, return closest time.
return_alts : bool
If `True`
Returns
-------
'''
if rise_set == 'rising':
condition = (alt[:-1] < horizon*unit) * (alt[1:] > horizon*unit)
elif rise_set == 'setting':
condition = (alt[:-1] > horizon*unit) * (alt[1:] < horizon*unit)
else:
raise ValueError('`rise_set` keyword may only be "rising" or '
'"setting"')
if sum(condition) < 1:
raise ValueError('Target does not rise/set within 24 hours')
if return_limits:
nearest_index = np.argwhere(condition)[0][0]
if alt[nearest_index] > horizon*unit:
lower_limit, upper_limit = t[nearest_index-1], t[nearest_index]
lower_alt, upper_alt = alt[nearest_index-1].value, alt[nearest_index].value
else:
lower_limit, upper_limit = t[nearest_index], t[nearest_index+1]
lower_alt, upper_alt = alt[nearest_index].value, alt[nearest_index+1].value
if return_alts:
return (lower_limit, upper_limit), (lower_alt, upper_alt)
else:
return lower_limit, upper_limit
return t[condition][0]
def _calc_altitude(time, target, location):
'''
Calculate altitude of `target` at `time` from `location`.
'''
return target.transform_to(AltAz(location=location,
obstime=time)).alt
def _two_point_interp(times, altitudes):
'''
Do linear interpolation between two `altitudes` at
two `times`.
Parameters
----------
times : `~astropy.time.Time`
altitudes : array of floats
Returns
-------
`~astropy.time.Time`
'''
slope = (altitudes[1] - altitudes[0])/(times[1].jd - times[0].jd)
return Time(times[1].jd - altitudes[1]/slope, format='jd')
def _calc_riseset(time, target, prev_next, rise_set, location,
horizon, N=150):
'''
Calculate the time at next rise/set of `target`.
Parameters
----------
time : `~astropy.time.Time`
Time of observation
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
location : `~astropy.coordinates.EarthLocation`
Location of observer
horizon : float
Number of degrees above/below actual horizon to use
for calculating rise/set times (i.e.,
-6 deg horizon = civil twilight, etc.)
N : int
Number of altitudes to compute when searching for
rise or set.
Returns
-------
ret1 : `~astropy.time.Time`
Time of rise/set
'''
if prev_next == 'next':
times = time + np.linspace(0, 1, N)*u.day
else:
times = time + np.linspace(-1, 0, N)*u.day
altitudes = _calc_altitude(times, target, location)
horizon_crossing_limits = _horiz_cross(times, altitudes, rise_set,
horizon,
return_limits=True,
return_alts=True)
return _two_point_interp(*horizon_crossing_limits)
def sunrise(time, location, which='nearest', horizon=0):
'''
Compute time of the next/previous/nearest sunrise, where
sunrise is defined as when the sun goes from altitudes
below `horizon` to above `horizon`.
Parameters
----------
time : `~astropy.time.Time`
Time of observation
location : `~astropy.coordinates.EarthLocation`
Location of observer
which : str - can be "next", "previous", or "nearest"; defaults to "nearest"
Choose which sunrise relative to the present `time` would you
like to calculate
horizon: float (optional); defaults to 0
Compute sunrise for horizon offset relative to actual horizon,
in degrees (i.e. horizon=-6 is would compute civil twilight)
Returns
-------
`~astropy.time.Time`
'''
if which == 'next' or which == 'nearest':
next_rise = _calc_riseset(test_time, get_sun(test_time), 'next',
'rising', location, horizon)
if which == 'next':
return next_rise
if which == 'previous' or which == 'nearest':
previous_rise = _calc_riseset(test_time, get_sun(test_time),
'previous', 'rising', location,
horizon)
if which == 'previous':
return previous_rise
if which == 'nearest':
if abs(time - previous_rise) < abs(time - next_rise):
return previous_rise
else:
return next_rise
def sunset(time, location, which='nearest', horizon=0):
'''
Compute time of the next/previous/nearest sunset, where
sunset is defined as when the sun goes from altitudes
above `horizon` to below `horizon`.
Parameters
----------
time : `~astropy.time.Time`
Time of observation
location : `~astropy.coordinates.EarthLocation`
Location of observer
which : str - can be "next", "previous", or "nearest"; defaults to "nearest"
Choose which sunset relative to the present `time` would you
like to calculate
horizon: float (optional); defaults to 0
Compute sunset for horizon offset relative to actual horizon,
in degrees (i.e. horizon=-6 is would compute civil twilight)
Returns
-------
`~astropy.time.Time`
'''
if which == 'next' or which == 'nearest':
next_set = _calc_riseset(test_time, get_sun(test_time), 'next',
'setting', location, horizon)
if which == 'next':
return next_set
if which == 'previous' or which == 'nearest':
previous_set = _calc_riseset(test_time, get_sun(test_time),
'previous', 'setting', location,
horizon)
if which == 'previous':
return previous_set
if which == 'nearest':
if abs(time - previous_set) < abs(time - next_set):
return previous_set
else:
return next_set
location = EarthLocation.from_geodetic(0, 65*u.degree, 0)
test_time = Time('2015-06-25 22:08:00')
next_rising = sunrise(test_time, location, which='next')
next_setting = sunset(test_time, location, which='next')
prev_rising = sunrise(test_time, location, which='previous')
prev_setting = sunset(test_time, location, which='previous')
print('\n'.join(map('{.iso}'.format,
[prev_setting, prev_rising,
next_setting, next_rising])))
# (should print True if sun is up at test_time)
print(prev_setting < prev_rising < test_time < next_setting < next_rising)
In [70]:
%timeit sunrise(test_time, location, which='next')
%timeit sunrise(test_time, location, which='nearest')
In [71]:
import ephem
obs = ephem.Observer()
obs.lat = '65:00:00'
obs.lon = '00:00:00'
obs.elevation = 0
obs.date = test_time.datetime
obs.horizon = 0
obs.pressure = 0
sun = ephem.Sun()
print(obs.previous_setting(sun, use_center='True'))
print(obs.previous_rising(sun, use_center='True'))
print(obs.next_setting(sun, use_center='True'))
print(obs.next_rising(sun, use_center='True'))
In [ ]: