# Source code for pvlib.solarposition

"""
Calculate the solar position using a variety of methods/packages.
"""

# Contributors:
# Rob Andrews (@Calama-Consulting), Calama Consulting, 2014
# Will Holmgren (@wholmgren), University of Arizona, 2014
# Tony Lorenzo (@alorenzo175), University of Arizona, 2015
# Cliff hansen (@cwhanse), Sandia National Laboratories, 2018

import os
import datetime as dt
try:
except ImportError:
try:
except ImportError:
pass

import numpy as np
import pandas as pd
import scipy.optimize as so
import warnings
import datetime

from pvlib import atmosphere
from pvlib.tools import datetime_to_djd, djd_to_datetime

[docs]def get_solarposition(time, latitude, longitude,
altitude=None, pressure=None,
method='nrel_numpy',
temperature=12, **kwargs):
"""
A convenience wrapper for the solar position calculators.

Parameters
----------
time : pandas.DatetimeIndex
Must be localized or UTC will be assumed.

latitude : float
Latitude in decimal degrees. Positive north of equator, negative
to south.

longitude : float
Longitude in decimal degrees. Positive east of prime meridian,
negative to west.

altitude : float, optional
If not specified, computed from pressure. Assumed to be 0 m
if pressure is not supplied.

pressure : float, optional
If not specified, computed from altitude. Assumed to be 101325 Pa
if altitude is not supplied.

method : string, default 'nrel_numpy'
'nrel_numpy' uses an implementation of the NREL SPA algorithm
described in [1] (default, recommended): :py:func:spa_python

'nrel_numba' uses an implementation of the NREL SPA algorithm
described in [1], but also compiles the code first:
:py:func:spa_python

'pyephem' uses the PyEphem package: :py:func:pyephem

'ephemeris' uses the pvlib ephemeris code: :py:func:ephemeris

'nrel_c' uses the NREL SPA C code [3]: :py:func:spa_c

temperature : float, default 12
Degrees C.

kwargs
Other keywords are passed to the solar position function
specified by the method argument.

References
----------
.. [1] I. Reda and A. Andreas, Solar position algorithm for solar radiation
applications. Solar Energy, vol. 76, no. 5, pp. 577-589, 2004.

.. [2] I. Reda and A. Andreas, Corrigendum to Solar position algorithm for
solar radiation applications. Solar Energy, vol. 81, no. 6, p. 838,
2007.

.. [3] NREL SPA code: https://midcdmz.nrel.gov/spa/
"""

if altitude is None and pressure is None:
altitude = 0.
pressure = 101325.
elif altitude is None:
altitude = atmosphere.pres2alt(pressure)
elif pressure is None:
pressure = atmosphere.alt2pres(altitude)

method = method.lower()
if isinstance(time, dt.datetime):
time = pd.DatetimeIndex([time, ])

if method == 'nrel_c':
ephem_df = spa_c(time, latitude, longitude, pressure, temperature,
**kwargs)
elif method == 'nrel_numba':
ephem_df = spa_python(time, latitude, longitude, altitude,
pressure, temperature,
how='numba', **kwargs)
elif method == 'nrel_numpy':
ephem_df = spa_python(time, latitude, longitude, altitude,
pressure, temperature,
how='numpy', **kwargs)
elif method == 'pyephem':
ephem_df = pyephem(time, latitude, longitude,
altitude=altitude,
pressure=pressure,
temperature=temperature, **kwargs)
elif method == 'ephemeris':
ephem_df = ephemeris(time, latitude, longitude, pressure, temperature,
**kwargs)
else:
raise ValueError('Invalid solar position method')

return ephem_df

[docs]def spa_c(time, latitude, longitude, pressure=101325, altitude=0,
temperature=12, delta_t=67.0,
raw_spa_output=False):
r"""
Calculate the solar position using the C implementation of the NREL
SPA code.

The source files for this code are located in './spa_c_files/', along with
a README file which describes how the C code is wrapped in Python.
and used in accordance with it's license.

This function is slower and no more accurate than :py:func:spa_python.

Parameters
----------
time : pandas.DatetimeIndex
Must be localized or UTC will be assumed.
latitude : float
Latitude in decimal degrees. Positive north of equator, negative
to south.
longitude : float
Longitude in decimal degrees. Positive east of prime meridian,
negative to west.
pressure : float, default 101325
Pressure in Pascals
altitude : float, default 0
Height above sea level. [m]
temperature : float, default 12
Temperature in C
delta_t : float, default 67.0
Difference between terrestrial time and UT1.
USNO has previous values and predictions [3]_.
raw_spa_output : bool, default False
If true, returns the raw SPA output.

Returns
-------
DataFrame
The DataFrame will have the following columns:
elevation,
azimuth,
zenith,
apparent_elevation,
apparent_zenith.

References
----------
.. [1] NREL SPA reference: https://midcdmz.nrel.gov/spa/

Note: The timezone field in the SPA C files is replaced with
time_zone to avoid a nameclash with the function __timezone that is
redefined by Python>=3.5. This issue is
Python bug 24643 <https://bugs.python.org/issue24643>_.

.. [2] Delta T: https://en.wikipedia.org/wiki/%CE%94T_(timekeeping)

.. [3] USNO delta T: https://maia.usno.navy.mil/products/deltaT

--------
pyephem, spa_python, ephemeris
"""

# Added by Rob Andrews (@Calama-Consulting), Calama Consulting, 2014
# Edited by Will Holmgren (@wholmgren), University of Arizona, 2014
# Edited by Tony Lorenzo (@alorenzo175), University of Arizona, 2015

try:
from pvlib.spa_c_files.spa_py import spa_calc
except ImportError:
raise ImportError('Could not import built-in SPA calculator. ' +
'You may need to recompile the SPA code.')

# if localized, convert to UTC. otherwise, assume UTC.
try:
time_utc = time.tz_convert('UTC')
except TypeError:
time_utc = time

spa_out = []

for date in time_utc:
spa_out.append(spa_calc(year=date.year,
month=date.month,
day=date.day,
hour=date.hour,
minute=date.minute,
second=date.second,
time_zone=0,  # date uses utc time
latitude=latitude,
longitude=longitude,
elevation=altitude,
pressure=pressure / 100,
temperature=temperature,
delta_t=delta_t
))

spa_df = pd.DataFrame(spa_out, index=time)

if raw_spa_output:
# rename "time_zone" from raw output from spa_c_files.spa_py.spa_calc()
# to "timezone" to match the API of pvlib.solarposition.spa_c()
return spa_df.rename(columns={'time_zone': 'timezone'})
else:
dfout = pd.DataFrame({'azimuth': spa_df['azimuth'],
'apparent_zenith': spa_df['zenith'],
'apparent_elevation': spa_df['e'],
'elevation': spa_df['e0'],
'zenith': 90 - spa_df['e0']})

return dfout

def _spa_python_import(how):
"""Compile spa.py appropriately"""

from pvlib import spa

# check to see if the spa module was compiled with numba
using_numba = spa.USE_NUMBA

if how == 'numpy' and using_numba:
# the spa module was compiled to numba code, so we need to
# reload the module without compiling
# the PVLIB_USE_NUMBA env variable is used to tell the module
# to not compile with numba
os.environ['PVLIB_USE_NUMBA'] = '0'
del os.environ['PVLIB_USE_NUMBA']
elif how == 'numba' and not using_numba:
# The spa module was not compiled to numba code, so set
# PVLIB_USE_NUMBA so it does compile to numba on reload.
os.environ['PVLIB_USE_NUMBA'] = '1'
del os.environ['PVLIB_USE_NUMBA']
elif how != 'numba' and how != 'numpy':
raise ValueError("how must be either 'numba' or 'numpy'")

return spa

def _datetime_to_unixtime(dtindex):
# convert a pandas datetime index to unixtime, making sure to handle
# different pandas units (ns, us, etc) and time zones correctly
if dtindex.tz is not None:
# epoch is 1970-01-01 00:00 UTC, but we need to match the input tz
# for compatibility with older pandas versions (e.g. v1.3.5)
epoch = pd.Timestamp("1970-01-01", tz="UTC").tz_convert(dtindex.tz)
else:
epoch = pd.Timestamp("1970-01-01")

return np.array((dtindex - epoch) / pd.Timedelta("1s"))

[docs]def spa_python(time, latitude, longitude,
altitude=0, pressure=101325, temperature=12, delta_t=67.0,
"""
Calculate the solar position using a python implementation of the
NREL SPA algorithm.

The details of the NREL SPA algorithm are described in [1]_.

If numba is installed, the functions can be compiled to
machine code and the function can be multithreaded.
Without numba, the function evaluates via numpy with
a slight performance hit.

Parameters
----------
time : pandas.DatetimeIndex
Must be localized or UTC will be assumed.
latitude : float
Latitude in decimal degrees. Positive north of equator, negative
to south.
longitude : float
Longitude in decimal degrees. Positive east of prime meridian,
negative to west.
altitude : float, default 0
Distance above sea level.
pressure : int or float, optional, default 101325
avg. yearly air pressure in Pascals.
temperature : int or float, optional, default 12
avg. yearly air temperature in degrees C.
delta_t : float, optional, default 67.0
Difference between terrestrial time and UT1.
If delta_t is None, uses spa.calculate_deltat
using time.year and time.month from pandas.DatetimeIndex.
For most simulations the default delta_t is sufficient.
*Note: delta_t = None will break code using nrel_numba,
this will be fixed in a future version.*
The USNO has historical and forecasted delta_t [3]_.
atmos_refrac : float, optional
The approximate atmospheric refraction (in degrees)
at sunrise and sunset.
how : str, optional, default 'numpy'
Options are 'numpy' or 'numba'. If numba >= 0.17.0
is installed, how='numba' will compile the spa functions
to machine code and run them multithreaded.
numthreads : int, optional, default 4
Number of threads to use if how == 'numba'.

Returns
-------
DataFrame
The DataFrame will have the following columns:
apparent_zenith (degrees),
zenith (degrees),
apparent_elevation (degrees),
elevation (degrees),
azimuth (degrees),
equation_of_time (minutes).

References
----------
.. [1] I. Reda and A. Andreas, Solar position algorithm for solar
radiation applications. Solar Energy, vol. 76, no. 5, pp. 577-589, 2004.

.. [2] I. Reda and A. Andreas, Corrigendum to Solar position algorithm for
solar radiation applications. Solar Energy, vol. 81, no. 6, p. 838,
2007.

.. [3] USNO delta T:
https://maia.usno.navy.mil/products/deltaT

--------
pyephem, spa_c, ephemeris
"""

# Added by Tony Lorenzo (@alorenzo175), University of Arizona, 2015

lat = latitude
lon = longitude
elev = altitude
pressure = pressure / 100  # pressure must be in millibars for calculation

atmos_refract = atmos_refract or 0.5667

if not isinstance(time, pd.DatetimeIndex):
try:
time = pd.DatetimeIndex(time)
except (TypeError, ValueError):
time = pd.DatetimeIndex([time, ])

unixtime = _datetime_to_unixtime(time)

spa = _spa_python_import(how)

delta_t = delta_t or spa.calculate_deltat(time.year, time.month)

app_zenith, zenith, app_elevation, elevation, azimuth, eot = \
spa.solar_position(unixtime, lat, lon, elev, pressure, temperature,

result = pd.DataFrame({'apparent_zenith': app_zenith, 'zenith': zenith,
'apparent_elevation': app_elevation,
'elevation': elevation, 'azimuth': azimuth,
'equation_of_time': eot},
index=time)

return result

[docs]def sun_rise_set_transit_spa(times, latitude, longitude, how='numpy',
"""
Calculate the sunrise, sunset, and sun transit times using the
NREL SPA algorithm.

The details of the NREL SPA algorithm are described in [1]_.

If numba is installed, the functions can be compiled to
machine code and the function can be multithreaded.
Without numba, the function evaluates via numpy with
a slight performance hit.

Parameters
----------
times : pandas.DatetimeIndex
Must be localized to the timezone for latitude and longitude.
latitude : float
Latitude in degrees, positive north of equator, negative to south
longitude : float
Longitude in degrees, positive east of prime meridian, negative to west
how : str, optional, default 'numpy'
Options are 'numpy' or 'numba'. If numba >= 0.17.0
is installed, how='numba' will compile the spa functions
to machine code and run them multithreaded.
delta_t : float, optional, default 67.0
Difference between terrestrial time and UT1.
If delta_t is None, uses spa.calculate_deltat
using times.year and times.month from pandas.DatetimeIndex.
For most simulations the default delta_t is sufficient.
*Note: delta_t = None will break code using nrel_numba,
this will be fixed in a future version.*
numthreads : int, optional, default 4
Number of threads to use if how == 'numba'.

Returns
-------
pandas.DataFrame
index is the same as input times argument
columns are 'sunrise', 'sunset', and 'transit'

References
----------
.. [1] Reda, I., Andreas, A., 2003. Solar position algorithm for solar
radiation applications. Technical report: NREL/TP-560- 34302. Golden,
USA, http://www.nrel.gov.
"""
# Added by Tony Lorenzo (@alorenzo175), University of Arizona, 2015

lat = latitude
lon = longitude

# times must be localized
if times.tz:
tzinfo = times.tz
else:
raise ValueError('times must be localized')

# must convert to midnight UTC on day of interest
utcday = pd.DatetimeIndex(times.date).tz_localize('UTC')
unixtime = _datetime_to_unixtime(utcday)

spa = _spa_python_import(how)

delta_t = delta_t or spa.calculate_deltat(times.year, times.month)

transit, sunrise, sunset = spa.transit_sunrise_sunset(

# arrays are in seconds since epoch format, need to conver to timestamps
transit = pd.to_datetime(transit*1e9, unit='ns', utc=True).tz_convert(
tzinfo).tolist()
sunrise = pd.to_datetime(sunrise*1e9, unit='ns', utc=True).tz_convert(
tzinfo).tolist()
sunset = pd.to_datetime(sunset*1e9, unit='ns', utc=True).tz_convert(
tzinfo).tolist()

return pd.DataFrame(index=times, data={'sunrise': sunrise,
'sunset': sunset,
'transit': transit})

def _ephem_convert_to_seconds_and_microseconds(date):
# utility from unreleased PyEphem 3.6.7.1
"""Converts a PyEphem date into seconds"""
microseconds = int(round(24 * 60 * 60 * 1000000 * date))
seconds, microseconds = divmod(microseconds, 1000000)
seconds -= 2209032000  # difference between epoch 1900 and epoch 1970
return seconds, microseconds

def _ephem_to_timezone(date, tzinfo):
# utility from unreleased PyEphem 3.6.7.1
""""Convert a PyEphem Date into a timezone aware python datetime"""
seconds, microseconds = _ephem_convert_to_seconds_and_microseconds(date)
date = dt.datetime.fromtimestamp(seconds, tzinfo)
date = date.replace(microsecond=microseconds)
return date

def _ephem_setup(latitude, longitude, altitude, pressure, temperature,
horizon):
import ephem
# initialize a PyEphem observer
obs = ephem.Observer()
obs.lat = str(latitude)
obs.lon = str(longitude)
obs.elevation = altitude
obs.pressure = pressure / 100.  # convert to mBar
obs.temp = temperature
obs.horizon = horizon

# the PyEphem sun
sun = ephem.Sun()
return obs, sun

[docs]def sun_rise_set_transit_ephem(times, latitude, longitude,
next_or_previous='next',
altitude=0,
pressure=101325,
temperature=12, horizon='0:00'):
"""
Calculate the next sunrise and sunset times using the PyEphem package.

Parameters
----------
time : pandas.DatetimeIndex
Must be localized
latitude : float
Latitude in degrees, positive north of equator, negative to south
longitude : float
Longitude in degrees, positive east of prime meridian, negative to west
next_or_previous : str
'next' or 'previous' sunrise and sunset relative to time
altitude : float, default 0
distance above sea level in meters.
pressure : int or float, optional, default 101325
air pressure in Pascals.
temperature : int or float, optional, default 12
air temperature in degrees C.
horizon : string, format +/-X:YY
arc degrees:arc minutes from geometrical horizon for sunrise and
sunset, e.g., horizon='+0:00' to use sun center crossing the
geometrical horizon to define sunrise and sunset,
horizon='-0:34' for when the sun's upper edge crosses the
geometrical horizon

Returns
-------
pandas.DataFrame
index is the same as input time argument
columns are 'sunrise', 'sunset', and 'transit'

--------
pyephem
"""

try:
import ephem
except ImportError:
raise ImportError('PyEphem must be installed')

# times must be localized
if times.tz:
tzinfo = times.tz
else:
raise ValueError('times must be localized')

obs, sun = _ephem_setup(latitude, longitude, altitude,
pressure, temperature, horizon)
# create lists of sunrise and sunset time localized to time.tz
if next_or_previous.lower() == 'next':
rising = obs.next_rising
setting = obs.next_setting
transit = obs.next_transit
elif next_or_previous.lower() == 'previous':
rising = obs.previous_rising
setting = obs.previous_setting
transit = obs.previous_transit
else:
raise ValueError("next_or_previous must be either 'next' or" +
" 'previous'")

sunrise = []
sunset = []
trans = []
for thetime in times:
thetime = thetime.to_pydatetime()
# older versions of pyephem ignore timezone when converting to its
# internal datetime format, so convert to UTC here to support
# all versions.  GH #1449
obs.date = ephem.Date(thetime.astimezone(datetime.timezone.utc))
sunrise.append(_ephem_to_timezone(rising(sun), tzinfo))
sunset.append(_ephem_to_timezone(setting(sun), tzinfo))
trans.append(_ephem_to_timezone(transit(sun), tzinfo))

return pd.DataFrame(index=times, data={'sunrise': sunrise,
'sunset': sunset,
'transit': trans})

[docs]def pyephem(time, latitude, longitude, altitude=0, pressure=101325,
temperature=12, horizon='+0:00'):
"""
Calculate the solar position using the PyEphem package.

Parameters
----------
time : pandas.DatetimeIndex
Must be localized or UTC will be assumed.
latitude : float
Latitude in decimal degrees. Positive north of equator, negative
to south.
longitude : float
Longitude in decimal degrees. Positive east of prime meridian,
negative to west.
altitude : float, default 0
Height above sea level in meters. [m]
pressure : int or float, optional, default 101325
air pressure in Pascals.
temperature : int or float, optional, default 12
air temperature in degrees C.
horizon : string, optional, default '+0:00'
arc degrees:arc minutes from geometrical horizon for sunrise and
sunset, e.g., horizon='+0:00' to use sun center crossing the
geometrical horizon to define sunrise and sunset,
horizon='-0:34' for when the sun's upper edge crosses the
geometrical horizon

Returns
-------
pandas.DataFrame
index is the same as input time argument
The DataFrame will have the following columns:
apparent_elevation, elevation,
apparent_azimuth, azimuth,
apparent_zenith, zenith.

--------
spa_python, spa_c, ephemeris
"""

# Written by Will Holmgren (@wholmgren), University of Arizona, 2014
try:
import ephem
except ImportError:
raise ImportError('PyEphem must be installed')

# if localized, convert to UTC. otherwise, assume UTC.
try:
time_utc = time.tz_convert('UTC')
except TypeError:
time_utc = time

sun_coords = pd.DataFrame(index=time)

obs, sun = _ephem_setup(latitude, longitude, altitude,
pressure, temperature, horizon)

# make and fill lists of the sun's altitude and azimuth
# this is the pressure and temperature corrected apparent alt/az.
alts = []
azis = []
for thetime in time_utc:
obs.date = ephem.Date(thetime)
sun.compute(obs)
alts.append(sun.alt)
azis.append(sun.az)

sun_coords['apparent_elevation'] = alts
sun_coords['apparent_azimuth'] = azis

# redo it for p=0 to get no atmosphere alt/az
obs.pressure = 0
alts = []
azis = []
for thetime in time_utc:
obs.date = ephem.Date(thetime)
sun.compute(obs)
alts.append(sun.alt)
azis.append(sun.az)

sun_coords['elevation'] = alts
sun_coords['azimuth'] = azis

# convert to degrees. add zenith
sun_coords['apparent_zenith'] = 90 - sun_coords['apparent_elevation']
sun_coords['zenith'] = 90 - sun_coords['elevation']

return sun_coords

[docs]def ephemeris(time, latitude, longitude, pressure=101325, temperature=12):
"""
Python-native solar position calculator.
The accuracy of this code is not guaranteed.
Consider using the built-in spa_c code or the PyEphem library.

Parameters
----------
time : pandas.DatetimeIndex
Must be localized or UTC will be assumed.
latitude : float
Latitude in decimal degrees. Positive north of equator, negative
to south.
longitude : float
Longitude in decimal degrees. Positive east of prime meridian,
negative to west.
pressure : float or Series, default 101325
Ambient pressure (Pascals)
temperature : float or Series, default 12
Ambient temperature (C)

Returns
-------

DataFrame with the following columns:

* apparent_elevation : apparent sun elevation accounting for
atmospheric refraction.
* elevation : actual elevation (not accounting for refraction)
of the sun in decimal degrees, 0 = on horizon.
The complement of the zenith angle.
* azimuth : Azimuth of the sun in decimal degrees East of North.
This is the complement of the apparent zenith angle.
* apparent_zenith : apparent sun zenith accounting for atmospheric
refraction.
* zenith : Solar zenith angle
* solar_time : Solar time in decimal hours (solar noon is 12.00).

References
-----------

.. [1] Grover Hughes' class and related class materials on Engineering
Astronomy at Sandia National Laboratories, 1985.

--------
pyephem, spa_c, spa_python

"""

# Added by Rob Andrews (@Calama-Consulting), Calama Consulting, 2014
# Edited by Will Holmgren (@wholmgren), University of Arizona, 2014

# Most comments in this function are from PVLIB_MATLAB or from
# pvlib-python's attempt to understand and fix problems with the
# algorithm. The comments are *not* based on the reference material.
# This helps a little bit:
# http://www.cv.nrao.edu/~rfisher/Ephemerides/times.html

# the inversion of longitude is due to the fact that this code was
# originally written for the convention that positive longitude were for
# locations west of the prime meridian. However, the correct convention (as
# of 2009) is to use negative longitudes for locations west of the prime
# meridian. Therefore, the user should input longitude values under the
# correct convention (e.g. Albuquerque is at -106 longitude), but it needs
# to be inverted for use in the code.

Latitude = latitude
Longitude = -1 * longitude

Abber = 20 / 3600.

# the SPA algorithm needs time to be expressed in terms of
# decimal UTC hours of the day of the year.

# if localized, convert to UTC. otherwise, assume UTC.
try:
time_utc = time.tz_convert('UTC')
except TypeError:
time_utc = time

# strip out the day of the year and calculate the decimal hour
DayOfYear = time_utc.dayofyear
DecHours = (time_utc.hour + time_utc.minute/60. + time_utc.second/3600. +
time_utc.microsecond/3600.e6)

# np.array needed for pandas > 0.20
UnivDate = np.array(DayOfYear)
UnivHr = np.array(DecHours)

Yr = np.array(time_utc.year) - 1900
YrBegin = 365 * Yr + np.floor((Yr - 1) / 4.) - 0.5

Ezero = YrBegin + UnivDate
T = Ezero / 36525.

# Calculate Greenwich Mean Sidereal Time (GMST)
GMST0 = 6 / 24. + 38 / 1440. + (
45.836 + 8640184.542 * T + 0.0929 * T ** 2) / 86400.
GMST0 = 360 * (GMST0 - np.floor(GMST0))
GMSTi = np.mod(GMST0 + 360 * (1.0027379093 * UnivHr / 24.), 360)

# Local apparent sidereal time
LocAST = np.mod((360 + GMSTi - Longitude), 360)

EpochDate = Ezero + UnivHr / 24.
T1 = EpochDate / 36525.

23.452294 - 0.0130125 * T1 - 1.64e-06 * T1 ** 2 + 5.03e-07 * T1 ** 3)
MlPerigee = 281.22083 + 4.70684e-05 * EpochDate + 0.000453 * T1 ** 2 + (
3e-06 * T1 ** 3)
MeanAnom = np.mod((358.47583 + 0.985600267 * EpochDate - 0.00015 *
T1 ** 2 - 3e-06 * T1 ** 3), 360)
Eccen = 0.01675104 - 4.18e-05 * T1 - 1.26e-07 * T1 ** 2
EccenAnom = MeanAnom
E = 0

while np.max(abs(EccenAnom - E)) > 0.0001:
E = EccenAnom

TrueAnom = (
2 * np.mod(np.degrees(np.arctan2(((1 + Eccen) / (1 - Eccen)) ** 0.5 *
EcLon = np.mod(MlPerigee + TrueAnom, 360) - Abber
DecR = np.arcsin(np.sin(ObliquityR)*np.sin(EcLonR))

RtAscen = np.degrees(np.arctan2(np.cos(ObliquityR)*np.sin(EcLonR),
np.cos(EcLonR)))

HrAngle = LocAST - RtAscen
HrAngle = HrAngle - (360 * (abs(HrAngle) > 180))

SunAz = np.degrees(np.arctan2(-np.sin(HrAngleR),
np.cos(LatR)*np.tan(DecR) -
np.sin(LatR)*np.cos(HrAngleR)))
SunAz[SunAz < 0] += 360

SunEl = np.degrees(np.arcsin(
np.cos(LatR) * np.cos(DecR) * np.cos(HrAngleR) +
np.sin(LatR) * np.sin(DecR)))

SolarTime = (180 + HrAngle) / 15.

# Calculate refraction correction
Elevation = SunEl
Refract = pd.Series(0., index=time_utc)

Refract[(Elevation > 5) & (Elevation <= 85)] = (
58.1/TanEl - 0.07/(TanEl**3) + 8.6e-05/(TanEl**5))

Refract[(Elevation > -0.575) & (Elevation <= 5)] = (
Elevation *
(-518.2 + Elevation*(103.4 + Elevation*(-12.79 + Elevation*0.711))) +
1735)

Refract[(Elevation > -1) & (Elevation <= -0.575)] = -20.774 / TanEl

Refract *= (283/(273. + temperature)) * (pressure/101325.) / 3600.

ApparentSunEl = SunEl + Refract

# make output DataFrame
DFOut = pd.DataFrame(index=time_utc)
DFOut['apparent_elevation'] = ApparentSunEl
DFOut['elevation'] = SunEl
DFOut['azimuth'] = SunAz
DFOut['apparent_zenith'] = 90 - ApparentSunEl
DFOut['zenith'] = 90 - SunEl
DFOut['solar_time'] = SolarTime
DFOut.index = time

return DFOut

[docs]def calc_time(lower_bound, upper_bound, latitude, longitude, attribute, value,
altitude=0, pressure=101325, temperature=12, horizon='+0:00',
xtol=1.0e-12):
"""
Calculate the time between lower_bound and upper_bound
where the attribute is equal to value. Uses PyEphem for
solar position calculations.

Parameters
----------
lower_bound : datetime.datetime
upper_bound : datetime.datetime
latitude : float
Latitude in decimal degrees. Positive north of equator, negative
to south.
longitude : float
Longitude in decimal degrees. Positive east of prime meridian,
negative to west.
attribute : str
The attribute of a pyephem.Sun object that
you want to solve for. Likely options are 'alt'
and 'az' (which must be given in radians).
value : int or float
The value of the attribute to solve for
altitude : float, default 0
Distance above sea level.
pressure : int or float, optional, default 101325
Air pressure in Pascals. Set to 0 for no
atmospheric correction.
temperature : int or float, optional, default 12
Air temperature in degrees C.
horizon : string, optional, default '+0:00'
arc degrees:arc minutes from geometrical horizon for sunrise and
sunset, e.g., horizon='+0:00' to use sun center crossing the
geometrical horizon to define sunrise and sunset,
horizon='-0:34' for when the sun's upper edge crosses the
geometrical horizon
xtol : float, optional, default 1.0e-12
The allowed error in the result from value

Returns
-------
datetime.datetime

Raises
------
ValueError
If the value is not contained between the bounds.
AttributeError
If the given attribute is not an attribute of a
PyEphem.Sun object.
"""
obs, sun = _ephem_setup(latitude, longitude, altitude,
pressure, temperature, horizon)

def compute_attr(thetime, target, attr):
obs.date = thetime
sun.compute(obs)
return getattr(sun, attr) - target

lb = datetime_to_djd(lower_bound)
ub = datetime_to_djd(upper_bound)

djd_root = so.brentq(compute_attr, lb, ub,
(value, attribute), xtol=xtol)

return djd_to_datetime(djd_root)

[docs]def pyephem_earthsun_distance(time):
"""
Calculates the distance from the earth to the sun using pyephem.

Parameters
----------
time : pandas.DatetimeIndex
Must be localized or UTC will be assumed.

Returns
-------
pd.Series. Earth-sun distance in AU.
"""

import ephem

sun = ephem.Sun()
earthsun = []
for thetime in time:
sun.compute(ephem.Date(thetime))
earthsun.append(sun.earth_distance)

return pd.Series(earthsun, index=time)

"""
Calculates the distance from the earth to the sun using the
NREL SPA algorithm.

The details of the NREL SPA algorithm are described in [1]_.

Parameters
----------
time : pandas.DatetimeIndex
Must be localized or UTC will be assumed.

how : str, optional, default 'numpy'
Options are 'numpy' or 'numba'. If numba >= 0.17.0
is installed, how='numba' will compile the spa functions
to machine code and run them multithreaded.

delta_t : float, optional, default 67.0
Difference between terrestrial time and UT1.
If delta_t is None, uses spa.calculate_deltat
using time.year and time.month from pandas.DatetimeIndex.
For most simulations the default delta_t is sufficient.
*Note: delta_t = None will break code using nrel_numba,
this will be fixed in a future version.*

numthreads : int, optional, default 4
Number of threads to use if how == 'numba'.

Returns
-------
dist : pd.Series
Earth-sun distance in AU.

References
----------
.. [1] Reda, I., Andreas, A., 2003. Solar position algorithm for solar
radiation applications. Technical report: NREL/TP-560- 34302. Golden,
USA, http://www.nrel.gov.
"""

if not isinstance(time, pd.DatetimeIndex):
try:
time = pd.DatetimeIndex(time)
except (TypeError, ValueError):
time = pd.DatetimeIndex([time, ])

unixtime = _datetime_to_unixtime(time)

spa = _spa_python_import(how)

delta_t = delta_t or spa.calculate_deltat(time.year, time.month)

dist = pd.Series(dist, index=time)

return dist

def _calculate_simple_day_angle(dayofyear, offset=1):
"""
Calculates the day angle for the Earth's orbit around the Sun.

Parameters
----------
dayofyear : numeric
offset : int, default 1
For the Spencer method, offset=1; for the ASCE method, offset=0

Returns
-------
day_angle : numeric
"""
return (2. * np.pi / 365.) * (dayofyear - offset)

[docs]def equation_of_time_spencer71(dayofyear):
"""
Equation of time from Duffie & Beckman and attributed to Spencer
(1971) and Iqbal (1983).

The coefficients correspond to the online copy of the Fourier
paper_ [1]_ in the Sundial Mailing list that was posted in 1998 by
Mac Oglesby from his correspondence with Macquarie University Prof.
John Pickard who added the following note.

was trying to use a hand calculator for calculating solar positions,
etc. He was extremely helpful and gave me a reprint of this paper. He
also pointed out an error in the original: in the series for E, the
constant was printed as 0.000075 rather than 0.0000075. I have
corrected the error in this version.

There appears to be another error in formula as printed in both
Duffie & Beckman's [2]_ and Frank Vignola's [3]_ books in which the
coefficient 0.04089 is printed instead of 0.040849, corresponding to
the value used in the Bird Clear Sky model implemented by Daryl
Myers [4]_ and printed in both the Fourier paper from the Sundial
Mailing List and R. Hulstrom's [5]_ book.

.. _Fourier paper: http://www.mail-archive.com/sundial@uni-koeln.de/msg01050.html

Parameters
----------
dayofyear : numeric

Returns
-------
equation_of_time : numeric
Difference in time between solar time and mean solar time in minutes.

References
----------
.. [1] J. W. Spencer, "Fourier series representation of the position of the
sun" in Search 2 (5), p. 172 (1971)

.. [2] J. A. Duffie and W. A. Beckman,  "Solar Engineering of Thermal
Processes, 3rd Edition" pp. 9-11, J. Wiley and Sons, New York (2006)

.. [3] Frank Vignola et al., "Solar And Infrared Radiation Measurements",
p. 13, CRC Press (2012)

.. [4] Daryl R. Myers, "Solar Radiation: Practical Modeling for Renewable
Energy Applications", p. 5 CRC Press (2013)

.. [5] Roland Hulstrom, "Solar Resources" p. 66, MIT Press (1989)

--------
equation_of_time_pvcdrom
"""
day_angle = _calculate_simple_day_angle(dayofyear)
# convert from radians to minutes per day = 24[h/day] * 60[min/h] / 2 / pi
eot = (1440.0 / 2 / np.pi) * (
0.0000075 +
0.001868 * np.cos(day_angle) - 0.032077 * np.sin(day_angle) -
0.014615 * np.cos(2.0 * day_angle) - 0.040849 * np.sin(2.0 * day_angle)
)
return eot

[docs]def equation_of_time_pvcdrom(dayofyear):
"""
Equation of time from PVCDROM.

PVCDROM_ is a website by Solar Power Lab at Arizona State
University (ASU)

.. _PVCDROM: http://www.pveducation.org/pvcdrom/2-properties-sunlight/solar-time

Parameters
----------
dayofyear : numeric

Returns
-------
equation_of_time : numeric
Difference in time between solar time and mean solar time in minutes.

References
----------
.. [1] Soteris A. Kalogirou, "Solar Energy Engineering Processes and
Systems, 2nd Edition" Elselvier/Academic Press (2009).

--------
equation_of_time_spencer71
"""
# day angle relative to Vernal Equinox, typically March 22 (day number 81)
bday = \
_calculate_simple_day_angle(dayofyear) - (2.0 * np.pi / 365.0) * 80.0
# same value but about 2x faster than Spencer (1971)
return 9.87 * np.sin(2.0 * bday) - 7.53 * np.cos(bday) - 1.5 * np.sin(bday)

[docs]def declination_spencer71(dayofyear):
"""
Solar declination from Duffie & Beckman and attributed to
Spencer (1971) and Iqbal (1983).

See [1]_ for details.

.. warning::
Return units are radians, not degrees.

Parameters
----------
dayofyear : numeric

Returns
-------
Angular position of the sun at solar noon relative to the plane of the
equator, approximately between +/-23.45 (degrees).

References
----------
.. [1] J. A. Duffie and W. A. Beckman,  "Solar Engineering of Thermal
Processes, 3rd Edition" pp. 13-14, J. Wiley and Sons, New York (2006)

.. [2] J. W. Spencer, "Fourier series representation of the position of the
sun" in Search 2 (5), p. 172 (1971)

.. [3] Daryl R. Myers, "Solar Radiation: Practical Modeling for Renewable
Energy Applications", p. 4 CRC Press (2013)

--------
declination_cooper69
"""
day_angle = _calculate_simple_day_angle(dayofyear)
return (
0.006918 -
0.399912 * np.cos(day_angle) + 0.070257 * np.sin(day_angle) -
0.006758 * np.cos(2. * day_angle) + 0.000907 * np.sin(2. * day_angle) -
0.002697 * np.cos(3. * day_angle) + 0.00148 * np.sin(3. * day_angle)
)

[docs]def declination_cooper69(dayofyear):
"""
Solar declination from Duffie & Beckman and attributed to Cooper (1969).

See [1]_ for details.

.. warning::
Return units are radians, not degrees.

Declination can be expressed using either sine or cosine:

.. math::

\\delta = 23.45 \\sin \\left( \\frac{2 \\pi}{365} \\left(n_{day} + 284
\\right) \\right) = -23.45 \\cos \\left( \\frac{2 \\pi}{365}
\\left(n_{day} + 10 \\right) \\right)

Parameters
----------
dayofyear : numeric

Returns
-------
Angular position of the sun at solar noon relative to the plane of the
equator, approximately between +/-23.45 (degrees).

References
----------
.. [1] J. A. Duffie and W. A. Beckman,  "Solar Engineering of Thermal
Processes, 3rd Edition" pp. 13-14, J. Wiley and Sons, New York (2006)

.. [2] J. H. Seinfeld and S. N. Pandis, "Atmospheric Chemistry and Physics"
p. 129, J. Wiley (1998)

.. [3] Daryl R. Myers, "Solar Radiation: Practical Modeling for Renewable
Energy Applications", p. 4 CRC Press (2013)

--------
declination_spencer71
"""
day_angle = _calculate_simple_day_angle(dayofyear)
dec = np.deg2rad(23.45 * np.sin(day_angle + (2.0 * np.pi / 365.0) * 285.0))
return dec

[docs]def solar_azimuth_analytical(latitude, hourangle, declination, zenith):
"""
Analytical expression of solar azimuth angle based on spherical
trigonometry.

Parameters
----------
latitude : numeric
hourangle : numeric
Hour angle in the local solar time in radians.
declination : numeric
Declination of the sun in radians.
zenith : numeric

Returns
-------
azimuth : numeric

References
----------
.. [1] J. A. Duffie and W. A. Beckman,  "Solar Engineering of Thermal
Processes, 3rd Edition" pp. 14, J. Wiley and Sons, New York (2006)

.. [2] J. H. Seinfeld and S. N. Pandis, "Atmospheric Chemistry and Physics"
p. 132, J. Wiley (1998)

.. [3] Wikipedia: Solar Azimuth Angle
<https://en.wikipedia.org/wiki/Solar_azimuth_angle>_

.. [4] PVCDROM: Azimuth Angle <http://www.pveducation.org/pvcdrom/2-
properties-sunlight/azimuth-angle>_

--------
declination_spencer71
declination_cooper69
hour_angle
solar_zenith_analytical
"""

numer = (np.cos(zenith) * np.sin(latitude) - np.sin(declination))
denom = (np.sin(zenith) * np.cos(latitude))

# cases that would generate new NaN values are safely ignored here
# since they are dealt with further below
with np.errstate(invalid='ignore', divide='ignore'):
cos_azi = numer / denom

# when zero division occurs, use the limit value of the analytical
# expression
cos_azi = \
np.where(np.isclose(denom,    0.0, rtol=0.0, atol=1e-8),  1.0, cos_azi)

# when too many round-ups in floating point math take cos_azi beyond
# 1.0, use 1.0
cos_azi = \
np.where(np.isclose(cos_azi,  1.0, rtol=0.0, atol=1e-8),  1.0, cos_azi)
cos_azi = \
np.where(np.isclose(cos_azi, -1.0, rtol=0.0, atol=1e-8), -1.0, cos_azi)

# when NaN values occur in input, ignore and pass to output
with np.errstate(invalid='ignore'):
sign_ha = np.sign(hourangle)

return sign_ha * np.arccos(cos_azi) + np.pi

[docs]def solar_zenith_analytical(latitude, hourangle, declination):
"""
Analytical expression of solar zenith angle based on spherical
trigonometry.

.. warning:: The analytic form neglects the effect of atmospheric
refraction.

Parameters
----------
latitude : numeric
hourangle : numeric
Hour angle in the local solar time in radians.
declination : numeric
Declination of the sun in radians.

Returns
-------
zenith : numeric

References
----------
.. [1] J. A. Duffie and W. A. Beckman,  "Solar Engineering of Thermal
Processes, 3rd Edition" pp. 14, J. Wiley and Sons, New York (2006)

.. [2] J. H. Seinfeld and S. N. Pandis, "Atmospheric Chemistry and
Physics" p. 132, J. Wiley (1998)

.. [3] Daryl R. Myers, "Solar Radiation: Practical Modeling for
Renewable Energy Applications", p. 5 CRC Press (2013)

.. [4] Wikipedia: Solar Zenith Angle
<https://en.wikipedia.org/wiki/Solar_zenith_angle>_

.. [5] PVCDROM: Elevation Angle
<https://www.pveducation.org/pvcdrom/properties-of-sunlight/
elevation-angle>_

--------
declination_spencer71
declination_cooper69
hour_angle
"""
return np.arccos(
np.cos(declination) * np.cos(latitude) * np.cos(hourangle) +
np.sin(declination) * np.sin(latitude)
)

[docs]def hour_angle(times, longitude, equation_of_time):
"""
Hour angle in local solar time. Zero at local solar noon.

Parameters
----------
times : :class:pandas.DatetimeIndex
Corresponding timestamps, must be localized to the timezone for the
longitude.
longitude : numeric
Longitude in degrees
equation_of_time : numeric
Equation of time in minutes.

Returns
-------
hour_angle : numeric
Hour angle in local solar time in degrees.

References
----------
.. [1] J. A. Duffie and W. A. Beckman,  "Solar Engineering of Thermal
Processes, 3rd Edition" pp. 13, J. Wiley and Sons, New York (2006)

.. [2] J. H. Seinfeld and S. N. Pandis, "Atmospheric Chemistry and Physics"
p. 132, J. Wiley (1998)

.. [3] Daryl R. Myers, "Solar Radiation: Practical Modeling for Renewable
Energy Applications", p. 5 CRC Press (2013)

--------
equation_of_time_spencer71
equation_of_time_pvcdrom
"""
# hours - timezone = (times - normalized_times) - (naive_times - times)
if times.tz is None:
times = times.tz_localize('utc')
tzs = np.array([ts.utcoffset().total_seconds() for ts in times]) / 3600

hrs_minus_tzs = (times - times.normalize()) / pd.Timedelta('1h') - tzs

# ensure array return instead of a version-dependent pandas <T>Index
return np.asarray(
15. * (hrs_minus_tzs - 12.) + longitude + equation_of_time / 4.)

def _hour_angle_to_hours(times, hourangle, longitude, equation_of_time):
"""converts hour angles in degrees to hours as a numpy array"""
if times.tz is None:
times = times.tz_localize('utc')
tzs = np.array([ts.utcoffset().total_seconds() for ts in times]) / 3600
hours = (hourangle - longitude - equation_of_time / 4.) / 15. + 12. + tzs
return np.asarray(hours)

def _local_times_from_hours_since_midnight(times, hours):
"""
converts hours since midnight from an array of floats to localized times
"""
tz_info = times.tz  # pytz timezone info
naive_times = times.tz_localize(None)  # naive but still localized
# normalize local, naive times to previous midnight and add the hours until
# sunrise, sunset, and transit
return pd.DatetimeIndex(
naive_times.normalize() + pd.to_timedelta(hours, unit='h'), tz=tz_info)

def _times_to_hours_after_local_midnight(times):
"""convert local pandas datetime indices to array of hours as floats"""
times = times.tz_localize(None)
hrs = (times - times.normalize()) / pd.Timedelta('1h')
return np.array(hrs)

[docs]def sun_rise_set_transit_geometric(times, latitude, longitude, declination,
equation_of_time):
"""
Geometric calculation of solar sunrise, sunset, and transit.

.. warning:: The geometric calculation assumes a circular earth orbit with
the sun as a point source at its center, and neglects the effect of
atmospheric refraction on zenith. The error depends on location and
time of year but is of order 10 minutes.

Parameters
----------
times : pandas.DatetimeIndex
Corresponding timestamps, must be localized to the timezone for the
latitude and longitude.
latitude : float
Latitude in degrees, positive north of equator, negative to south
longitude : float
Longitude in degrees, positive east of prime meridian, negative to west
declination : numeric
declination angle in radians at times
equation_of_time : numeric
difference in time between solar time and mean solar time in minutes

Returns
-------
sunrise : datetime
localized sunrise time
sunset : datetime
localized sunset time
transit : datetime
localized sun transit time

References
----------
.. [1] J. A. Duffie and W. A. Beckman,  "Solar Engineering of Thermal
Processes, 3rd Edition," J. Wiley and Sons, New York (2006)

.. [2] Frank Vignola et al., "Solar And Infrared Radiation Measurements,"
CRC Press (2012)

"""