From da453db0551019647fe029b15fbfe0ed83660789 Mon Sep 17 00:00:00 2001 From: Tony Lorenzo Date: Mon, 23 Mar 2015 09:18:54 -0700 Subject: [PATCH 1/2] implement NREL SPA code in numpy or compiled numba --- .travis.yml | 1 + docs/sphinx/source/whatsnew/v0.1.0.txt | 1 + pvlib/__init__.py | 1 + pvlib/solarposition.py | 167 +++- pvlib/spa.py | 1081 ++++++++++++++++++++++++ pvlib/spa_c_files/cspa_py.pxd | 2 + pvlib/spa_c_files/spa_py.c | 148 ++-- pvlib/spa_c_files/spa_py.pyx | 8 +- pvlib/test/test_solarposition.py | 85 +- pvlib/test/test_spa.py | 279 ++++++ 10 files changed, 1694 insertions(+), 79 deletions(-) create mode 100644 pvlib/spa.py create mode 100644 pvlib/test/test_spa.py diff --git a/.travis.yml b/.travis.yml index 05e5b00e29..6940d752f9 100644 --- a/.travis.yml +++ b/.travis.yml @@ -26,6 +26,7 @@ before_install: install: - echo "install" - conda install --yes python=$TRAVIS_PYTHON_VERSION numpy scipy pandas=$PD_VERSION nose pytz ephem + - test $PD_VERSION != 0.13.1 && conda install --yes "numba>=0.17" || echo "Not installing numba" - conda install --yes coverage - pip install coveralls - python setup.py install diff --git a/docs/sphinx/source/whatsnew/v0.1.0.txt b/docs/sphinx/source/whatsnew/v0.1.0.txt index 0402aef0ea..e1da046bb6 100644 --- a/docs/sphinx/source/whatsnew/v0.1.0.txt +++ b/docs/sphinx/source/whatsnew/v0.1.0.txt @@ -24,6 +24,7 @@ New features * Library is Python 3.3 and 3.4 compatible * Add What's New section to docs (:issue:`10`) * Add PyEphem option to solar position calculations. +* Add a Python translation of NREL's SPA algorithm. * ``irradiance.py`` has more AOI, projection, and irradiance sum and calculation functions * TMY data import has a ``coerce_year`` option * TMY data can be loaded from a url (:issue:`5`) diff --git a/pvlib/__init__.py b/pvlib/__init__.py index 634f894160..20ed1d4d26 100644 --- a/pvlib/__init__.py +++ b/pvlib/__init__.py @@ -10,3 +10,4 @@ from pvlib import tmy from pvlib import tracking from pvlib import pvsystem +from pvlib import spa diff --git a/pvlib/solarposition.py b/pvlib/solarposition.py index 481da6fd50..24cf9c98e9 100644 --- a/pvlib/solarposition.py +++ b/pvlib/solarposition.py @@ -5,11 +5,20 @@ # Contributors: # Rob Andrews (@Calama-Consulting), Calama Consulting, 2014 # Will Holmgren (@wholmgren), University of Arizona, 2014 +# Tony Lorenzo (@alorenzo175), University of Arizona, 2015 from __future__ import division +import os import logging pvl_logger = logging.getLogger('pvlib') import datetime as dt +try: + from importlib import reload +except ImportError: + try: + from imp import reload + except ImportError: + pass import numpy as np @@ -32,7 +41,11 @@ def get_solarposition(time, location, method='pyephem', pressure=101325, method : string 'pyephem' uses the PyEphem package (default): :func:`pyephem` - 'spa' uses the spa code: :func:`spa` + 'spa' or 'spa_c' uses the spa code: :func:`spa` + + 'spa_numpy' uses SPA code translated to python: :func: `spa_python` + + 'spa_numba' uses SPA code translated to python: :func: `spa_python` 'ephemeris' uses the pvlib ephemeris code: :func:`ephemeris` pressure : float @@ -45,10 +58,16 @@ def get_solarposition(time, location, method='pyephem', pressure=101325, if isinstance(time, dt.datetime): time = pd.DatetimeIndex([time, ]) - if method == 'spa': - ephem_df = spa(time, location) + if method == 'spa' or method == 'spa_c': + ephem_df = spa(time, location, pressure, temperature) elif method == 'pyephem': ephem_df = pyephem(time, location, pressure, temperature) + elif method == 'spa_numpy': + ephem_df = spa_python(time, location, pressure, temperature, + how='numpy') + elif method == 'spa_numba': + ephem_df = spa_python(time, location, pressure, temperature, + how='numba') elif method == 'ephemeris': ephem_df = ephemeris(time, location, pressure, temperature) else: @@ -57,7 +76,8 @@ def get_solarposition(time, location, method='pyephem', pressure=101325, return ephem_df -def spa(time, location, raw_spa_output=False): +def spa(time, location, pressure=101325, temperature=12, delta_t=67.0, + raw_spa_output=False): ''' Calculate the solar position using the C implementation of the NREL SPA code @@ -69,6 +89,13 @@ def spa(time, location, raw_spa_output=False): ---------- time : pandas.DatetimeIndex location : pvlib.Location object + pressure : float + Pressure in Pascals + temperature : float + Temperature in C + delta_t : float + Difference between terrestrial time and UT1. + USNO has previous values and predictions. raw_spa_output : bool If true, returns the raw SPA output. @@ -78,15 +105,19 @@ def spa(time, location, raw_spa_output=False): The DataFrame will have the following columns: elevation, azimuth, - zenith. + zenith, + apparent_elevation, + apparent_zenith. References ---------- NREL SPA code: http://rredc.nrel.gov/solar/codesandalgorithms/spa/ + USNO delta T: http://www.usno.navy.mil/USNO/earth-orientation/eo-products/long-term ''' # 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 @@ -102,27 +133,129 @@ def spa(time, location, raw_spa_output=False): 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, - timezone=0, # timezone corrections handled above - latitude=location.latitude, - longitude=location.longitude, - elevation=location.altitude)) + month=date.month, + day=date.day, + hour=date.hour, + minute=date.minute, + second=date.second, + timezone=0, # timezone corrections handled above + latitude=location.latitude, + longitude=location.longitude, + elevation=location.altitude, + pressure=pressure / 100, + temperature=temperature, + delta_t=delta_t + )) spa_df = pd.DataFrame(spa_out, index=time_utc).tz_convert(location.tz) if raw_spa_output: return spa_df else: - dfout = spa_df[['zenith', 'azimuth']] - dfout['elevation'] = 90 - dfout.zenith - + 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(time, location, pressure=101325, temperature=12, delta_t=None, + atmos_refract=None, how='numpy', numthreads=4): + """ + Calculate the solar position using a python translation of the + NREL SPA code. + + 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 + location : pvlib.Location object + pressure : int or float, optional + avg. yearly air pressure in Pascals. + temperature : int or float, optional + avg. yearly air temperature in degrees C. + delta_t : float, optional + Difference between terrestrial time and UT1. + By default, use USNO historical data and predictions + how : str, optional + If numba >= 0.17.0 is installed, how='numba' will + compile the spa functions to machine code and + run them multithreaded. Otherwise, numpy calculations + are used. + numthreads : int, optional + Number of threads to use if how == 'numba'. + + Returns + ------- + DataFrame + The DataFrame will have the following columns: + apparent_zenith, zenith, + apparent_elevation, elevation, + azimuth. + + + 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. + """ + + # Added by Tony Lorenzo (@alorenzo175), University of Arizona, 2015 + + from pvlib import spa + pvl_logger.debug('Calculating solar position with spa_python code') + + lat = location.latitude + lon = location.longitude + elev = location.altitude + pressure = pressure / 100 # pressure must be in millibars for calculation + delta_t = delta_t or 67.0 + 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 = localize_to_utc(time, location).astype(int)/10**9 + + using_numba = spa.USE_NUMBA + if how == 'numpy' and using_numba: + os.environ['PVLIB_USE_NUMBA'] = '0' + pvl_logger.debug('Reloading spa module without compiling') + spa = reload(spa) + del os.environ['PVLIB_USE_NUMBA'] + elif how == 'numba' and not using_numba: + os.environ['PVLIB_USE_NUMBA'] = '1' + pvl_logger.debug('Reloading spa module, compiling with numba') + spa = reload(spa) + del os.environ['PVLIB_USE_NUMBA'] + elif how != 'numba' and how != 'numpy': + raise ValueError("how must be either 'numba' or 'numpy'") + + app_zenith, zenith, app_elevation, elevation, azimuth = spa.solar_position( + unixtime, lat, lon, elev, pressure, temperature, delta_t, atmos_refract, + numthreads) + result = pd.DataFrame({'apparent_zenith':app_zenith, 'zenith':zenith, + 'apparent_elevation':app_elevation, + 'elevation':elevation, 'azimuth':azimuth}, + index = time) + + try: + result = result.tz_convert(location.tz) + except TypeError: + result = result.tz_localize(location.tz) + + return result + + def _ephem_setup(location, pressure, temperature): import ephem # initialize a PyEphem observer diff --git a/pvlib/spa.py b/pvlib/spa.py new file mode 100644 index 0000000000..0dac087588 --- /dev/null +++ b/pvlib/spa.py @@ -0,0 +1,1081 @@ +""" +Calculate the solar position using a translation of NREL SPA code either using +numpy arrays or compiling the code to machine language with numba. +""" + +# Contributors: +# Created by Tony Lorenzo (@alorenzo175), Univ. of Arizona, 2015 + +from __future__ import division +import os +import threading +import warnings +import logging +pvl_logger = logging.getLogger('pvlib') + + +import numpy as np + + +def nocompile(*args, **kwargs): + return lambda func: func + + +if os.getenv('PVLIB_USE_NUMBA', '0') != '0': + try: + from numba import vectorize, jit, float64, guvectorize, __version__ + major, minor = __version__.split('.')[:2] + if int(major + minor) >= 17: + # need at least numba >= 0.17.0 + jcompile = jit + USE_NUMBA = True + else: + warnings.warn('Numba version must be >= 0.17.0, falling back to numpy') + jcompile = nocompile + USE_NUMBA = False + except ImportError: + warnings.warn('Could not import numba, falling back to numpy calculation') + jcompile = nocompile + USE_NUMBA = False +else: + jcompile = nocompile + USE_NUMBA = False + + +TABLE_1_DICT = { + 'L0':np.array( + [[175347046.0, 0.0, 0.0], + [3341656.0, 4.6692568, 6283.07585], + [34894.0, 4.6261, 12566.1517], + [3497.0, 2.7441, 5753.3849], + [3418.0, 2.8289, 3.5231], + [3136.0, 3.6277, 77713.7715], + [2676.0, 4.4181, 7860.4194], + [2343.0, 6.1352, 3930.2097], + [1324.0, 0.7425, 11506.7698], + [1273.0, 2.0371, 529.691], + [1199.0, 1.1096, 1577.3435], + [990.0, 5.233, 5884.927], + [902.0, 2.045, 26.298], + [857.0, 3.508, 398.149], + [780.0, 1.179, 5223.694], + [753.0, 2.533, 5507.553], + [505.0, 4.583, 18849.228], + [492.0, 4.205, 775.523], + [357.0, 2.92, 0.067], + [317.0, 5.849, 11790.629], + [284.0, 1.899, 796.298], + [271.0, 0.315, 10977.079], + [243.0, 0.345, 5486.778], + [206.0, 4.806, 2544.314], + [205.0, 1.869, 5573.143], + [202.0, 2.458, 6069.777], + [156.0, 0.833, 213.299], + [132.0, 3.411, 2942.463], + [126.0, 1.083, 20.775], + [115.0, 0.645, 0.98], + [103.0, 0.636, 4694.003], + [102.0, 0.976, 15720.839], + [102.0, 4.267, 7.114], + [99.0, 6.21, 2146.17], + [98.0, 0.68, 155.42], + [86.0, 5.98, 161000.69], + [85.0, 1.3, 6275.96], + [85.0, 3.67, 71430.7], + [80.0, 1.81, 17260.15], + [79.0, 3.04, 12036.46], + [75.0, 1.76, 5088.63], + [74.0, 3.5, 3154.69], + [74.0, 4.68, 801.82], + [70.0, 0.83, 9437.76], + [62.0, 3.98, 8827.39], + [61.0, 1.82, 7084.9], + [57.0, 2.78, 6286.6], + [56.0, 4.39, 14143.5], + [56.0, 3.47, 6279.55], + [52.0, 0.19, 12139.55], + [52.0, 1.33, 1748.02], + [51.0, 0.28, 5856.48], + [49.0, 0.49, 1194.45], + [41.0, 5.37, 8429.24], + [41.0, 2.4, 19651.05], + [39.0, 6.17, 10447.39], + [37.0, 6.04, 10213.29], + [37.0, 2.57, 1059.38], + [36.0, 1.71, 2352.87], + [36.0, 1.78, 6812.77], + [33.0, 0.59, 17789.85], + [30.0, 0.44, 83996.85], + [30.0, 2.74, 1349.87], + [25.0, 3.16, 4690.48]]), + 'L1':np.array( + [[628331966747.0, 0.0, 0.0], + [206059.0, 2.678235, 6283.07585], + [4303.0, 2.6351, 12566.1517], + [425.0, 1.59, 3.523], + [119.0, 5.796, 26.298], + [109.0, 2.966, 1577.344], + [93.0, 2.59, 18849.23], + [72.0, 1.14, 529.69], + [68.0, 1.87, 398.15], + [67.0, 4.41, 5507.55], + [59.0, 2.89, 5223.69], + [56.0, 2.17, 155.42], + [45.0, 0.4, 796.3], + [36.0, 0.47, 775.52], + [29.0, 2.65, 7.11], + [21.0, 5.34, 0.98], + [19.0, 1.85, 5486.78], + [19.0, 4.97, 213.3], + [17.0, 2.99, 6275.96], + [16.0, 0.03, 2544.31], + [16.0, 1.43, 2146.17], + [15.0, 1.21, 10977.08], + [12.0, 2.83, 1748.02], + [12.0, 3.26, 5088.63], + [12.0, 5.27, 1194.45], + [12.0, 2.08, 4694.0], + [11.0, 0.77, 553.57], + [10.0, 1.3, 6286.6], + [10.0, 4.24, 1349.87], + [9.0, 2.7, 242.73], + [9.0, 5.64, 951.72], + [8.0, 5.3, 2352.87], + [6.0, 2.65, 9437.76], + [6.0, 4.67, 4690.48]]), + 'L2':np.array( + [[52919.0, 0.0, 0.0], + [8720.0, 1.0721, 6283.0758], + [309.0, 0.867, 12566.152], + [27.0, 0.05, 3.52], + [16.0, 5.19, 26.3], + [16.0, 3.68, 155.42], + [10.0, 0.76, 18849.23], + [9.0, 2.06, 77713.77], + [7.0, 0.83, 775.52], + [5.0, 4.66, 1577.34], + [4.0, 1.03, 7.11], + [4.0, 3.44, 5573.14], + [3.0, 5.14, 796.3], + [3.0, 6.05, 5507.55], + [3.0, 1.19, 242.73], + [3.0, 6.12, 529.69], + [3.0, 0.31, 398.15], + [3.0, 2.28, 553.57], + [2.0, 4.38, 5223.69], + [2.0, 3.75, 0.98]]), + 'L3':np.array( + [[289.0, 5.844, 6283.076], + [35.0, 0.0, 0.0], + [17.0, 5.49, 12566.15], + [3.0, 5.2, 155.42], + [1.0, 4.72, 3.52], + [1.0, 5.3, 18849.23], + [1.0, 5.97, 242.73]]), + 'L4':np.array( + [[114.0, 3.142, 0.0], + [8.0, 4.13, 6283.08], + [1.0, 3.84, 12566.15]]), + 'L5':np.array( + [[1.0, 3.14, 0.0]]), + 'B0':np.array( + [[280.0, 3.199, 84334.662], + [102.0, 5.422, 5507.553], + [80.0, 3.88, 5223.69], + [44.0, 3.7, 2352.87], + [32.0, 4.0, 1577.34]]), + 'B1':np.array( + [[9.0, 3.9, 5507.55], + [6.0, 1.73, 5223.69]]), + 'R0':np.array( + [[100013989.0, 0.0, 0.0], + [1670700.0, 3.0984635, 6283.07585], + [13956.0, 3.05525, 12566.1517], + [3084.0, 5.1985, 77713.7715], + [1628.0, 1.1739, 5753.3849], + [1576.0, 2.8469, 7860.4194], + [925.0, 5.453, 11506.77], + [542.0, 4.564, 3930.21], + [472.0, 3.661, 5884.927], + [346.0, 0.964, 5507.553], + [329.0, 5.9, 5223.694], + [307.0, 0.299, 5573.143], + [243.0, 4.273, 11790.629], + [212.0, 5.847, 1577.344], + [186.0, 5.022, 10977.079], + [175.0, 3.012, 18849.228], + [110.0, 5.055, 5486.778], + [98.0, 0.89, 6069.78], + [86.0, 5.69, 15720.84], + [86.0, 1.27, 161000.69], + [65.0, 0.27, 17260.15], + [63.0, 0.92, 529.69], + [57.0, 2.01, 83996.85], + [56.0, 5.24, 71430.7], + [49.0, 3.25, 2544.31], + [47.0, 2.58, 775.52], + [45.0, 5.54, 9437.76], + [43.0, 6.01, 6275.96], + [39.0, 5.36, 4694.0], + [38.0, 2.39, 8827.39], + [37.0, 0.83, 19651.05], + [37.0, 4.9, 12139.55], + [36.0, 1.67, 12036.46], + [35.0, 1.84, 2942.46], + [33.0, 0.24, 7084.9], + [32.0, 0.18, 5088.63], + [32.0, 1.78, 398.15], + [28.0, 1.21, 6286.6], + [28.0, 1.9, 6279.55], + [26.0, 4.59, 10447.39]]), + 'R1':np.array( + [[103019.0, 1.10749, 6283.07585], + [1721.0, 1.0644, 12566.1517], + [702.0, 3.142, 0.0], + [32.0, 1.02, 18849.23], + [31.0, 2.84, 5507.55], + [25.0, 1.32, 5223.69], + [18.0, 1.42, 1577.34], + [10.0, 5.91, 10977.08], + [9.0, 1.42, 6275.96], + [9.0, 0.27, 5486.78]]), + 'R2':np.array( + [[4359.0, 5.7846, 6283.0758], + [124.0, 5.579, 12566.152], + [12.0, 3.14, 0.0], + [9.0, 3.63, 77713.77], + [6.0, 1.87, 5573.14], + [3.0, 5.47, 18849.23]]), + 'R3':np.array( + [[145.0, 4.273, 6283.076], + [7.0, 3.92, 12566.15]]), + 'R4':np.array( + [[4.0, 2.56, 6283.08]]) +} + + +TABLE_1_DICT['L1'].resize((64,3)) +TABLE_1_DICT['L2'].resize((64,3)) +TABLE_1_DICT['L3'].resize((64,3)) +TABLE_1_DICT['L4'].resize((64,3)) +TABLE_1_DICT['L5'].resize((64,3)) + +TABLE_1_DICT['B1'].resize((5,3)) + +TABLE_1_DICT['R1'].resize((40,3)) +TABLE_1_DICT['R2'].resize((40,3)) +TABLE_1_DICT['R3'].resize((40,3)) +TABLE_1_DICT['R4'].resize((40,3)) + + +HELIO_LONG_TABLE = np.array([TABLE_1_DICT['L0'], + TABLE_1_DICT['L1'], + TABLE_1_DICT['L2'], + TABLE_1_DICT['L3'], + TABLE_1_DICT['L4'], + TABLE_1_DICT['L5']]) + + +HELIO_LAT_TABLE = np.array([TABLE_1_DICT['B0'], + TABLE_1_DICT['B1']]) + + +HELIO_RADIUS_TABLE = np.array([TABLE_1_DICT['R0'], + TABLE_1_DICT['R1'], + TABLE_1_DICT['R2'], + TABLE_1_DICT['R3'], + TABLE_1_DICT['R4']]) + + +NUTATION_ABCD_ARRAY = np.array([ + [-171996,-174.2,92025,8.9], + [-13187,-1.6,5736,-3.1], + [-2274,-0.2,977,-0.5], + [2062,0.2,-895,0.5], + [1426,-3.4,54,-0.1], + [712,0.1,-7,0], + [-517,1.2,224,-0.6], + [-386,-0.4,200,0], + [-301,0,129,-0.1], + [217,-0.5,-95,0.3], + [-158,0,0,0], + [129,0.1,-70,0], + [123,0,-53,0], + [63,0,0,0], + [63,0.1,-33,0], + [-59,0,26,0], + [-58,-0.1,32,0], + [-51,0,27,0], + [48,0,0,0], + [46,0,-24,0], + [-38,0,16,0], + [-31,0,13,0], + [29,0,0,0], + [29,0,-12,0], + [26,0,0,0], + [-22,0,0,0], + [21,0,-10,0], + [17,-0.1,0,0], + [16,0,-8,0], + [-16,0.1,7,0], + [-15,0,9,0], + [-13,0,7,0], + [-12,0,6,0], + [11,0,0,0], + [-10,0,5,0], + [-8,0,3,0], + [7,0,-3,0], + [-7,0,0,0], + [-7,0,3,0], + [-7,0,3,0], + [6,0,0,0], + [6,0,-3,0], + [6,0,-3,0], + [-6,0,3,0], + [-6,0,3,0], + [5,0,0,0], + [-5,0,3,0], + [-5,0,3,0], + [-5,0,3,0], + [4,0,0,0], + [4,0,0,0], + [4,0,0,0], + [-4,0,0,0], + [-4,0,0,0], + [-4,0,0,0], + [3,0,0,0], + [-3,0,0,0], + [-3,0,0,0], + [-3,0,0,0], + [-3,0,0,0], + [-3,0,0,0], + [-3,0,0,0], + [-3,0,0,0], +]) + + +NUTATION_YTERM_ARRAY = np.array([ + [0,0,0,0,1], + [-2,0,0,2,2], + [0,0,0,2,2], + [0,0,0,0,2], + [0,1,0,0,0], + [0,0,1,0,0], + [-2,1,0,2,2], + [0,0,0,2,1], + [0,0,1,2,2], + [-2,-1,0,2,2], + [-2,0,1,0,0], + [-2,0,0,2,1], + [0,0,-1,2,2], + [2,0,0,0,0], + [0,0,1,0,1], + [2,0,-1,2,2], + [0,0,-1,0,1], + [0,0,1,2,1], + [-2,0,2,0,0], + [0,0,-2,2,1], + [2,0,0,2,2], + [0,0,2,2,2], + [0,0,2,0,0], + [-2,0,1,2,2], + [0,0,0,2,0], + [-2,0,0,2,0], + [0,0,-1,2,1], + [0,2,0,0,0], + [2,0,-1,0,1], + [-2,2,0,2,2], + [0,1,0,0,1], + [-2,0,1,0,1], + [0,-1,0,0,1], + [0,0,2,-2,0], + [2,0,-1,2,1], + [2,0,1,2,2], + [0,1,0,2,2], + [-2,1,1,0,0], + [0,-1,0,2,2], + [2,0,0,2,1], + [2,0,1,0,0], + [-2,0,2,2,2], + [-2,0,1,2,1], + [2,0,-2,0,1], + [2,0,0,0,1], + [0,-1,1,0,0], + [-2,-1,0,2,1], + [-2,0,0,0,1], + [0,0,2,2,1], + [-2,0,2,0,1], + [-2,1,0,2,1], + [0,0,1,-2,0], + [-1,0,1,0,0], + [-2,1,0,0,0], + [1,0,0,0,0], + [0,0,1,2,0], + [0,0,-2,2,2], + [-1,-1,1,0,0], + [0,1,1,0,0], + [0,-1,1,2,2], + [2,-1,-1,2,2], + [0,0,3,2,2], + [2,-1,0,2,2], +]) + + +@jcompile('float64(int64, int64, int64, int64, int64, int64, int64)', + nopython=True) +def julian_day_dt(year, month, day, hour, minute, second, microsecond): + """This is the original way to calculate the julian day from the NREL paper. + However, it is much faster to convert to unix/epoch time and then convert + to julian day. Note that the date must be UTC.""" + if month <= 2: + year = year-1 + month = month+12 + a = int(year/100) + b = 2 - a + int(a * 0.25) + frac_of_day = (microsecond + (second + minute * 60 + hour * 3600 ) + ) * 1.0 / (3600*24) + d = day + frac_of_day + jd = (int(365.25 * (year + 4716)) + int(30.6001 * (month + 1)) + d + + b - 1524.5) + return jd + + +@jcompile('float64(float64)', nopython=True) +def julian_day(unixtime): + jd = unixtime*1.0/86400 + 2440587.5 + return jd + + +@jcompile('float64(float64, float64)', nopython=True) +def julian_ephemeris_day(julian_day, delta_t): + jde = julian_day + delta_t * 1.0 / 86400 + return jde + + +@jcompile('float64(float64)', nopython=True) +def julian_century(julian_day): + jc = (julian_day - 2451545) * 1.0 / 36525 + return jc + + +@jcompile('float64(float64)', nopython=True) +def julian_ephemeris_century(julian_ephemeris_day): + jce = (julian_ephemeris_day - 2451545) * 1.0 / 36525 + return jce + + +@jcompile('float64(float64)', nopython=True) +def julian_ephemeris_millennium(julian_ephemeris_century): + jme = julian_ephemeris_century * 1.0 / 10 + return jme + + +@jcompile('float64(float64)', nopython=True) +def heliocentric_longitude(jme): + l0 = 0.0 + l1 = 0.0 + l2 = 0.0 + l3 = 0.0 + l4 = 0.0 + l5 = 0.0 + + for row in range(HELIO_LONG_TABLE.shape[1]): + l0 += (HELIO_LONG_TABLE[0, row, 0] + * np.cos(HELIO_LONG_TABLE[0, row, 1] + + HELIO_LONG_TABLE[0, row, 2] * jme) + ) + l1 += (HELIO_LONG_TABLE[1, row, 0] + * np.cos(HELIO_LONG_TABLE[1, row, 1] + + HELIO_LONG_TABLE[1, row, 2] * jme) + ) + l2 += (HELIO_LONG_TABLE[2, row, 0] + * np.cos(HELIO_LONG_TABLE[2, row, 1] + + HELIO_LONG_TABLE[2, row, 2] * jme) + ) + l3 += (HELIO_LONG_TABLE[3, row, 0] + * np.cos(HELIO_LONG_TABLE[3, row, 1] + + HELIO_LONG_TABLE[3, row, 2] * jme) + ) + l4 += (HELIO_LONG_TABLE[4, row, 0] + * np.cos(HELIO_LONG_TABLE[4, row, 1] + + HELIO_LONG_TABLE[4, row, 2] * jme) + ) + l5 += (HELIO_LONG_TABLE[5, row, 0] + * np.cos(HELIO_LONG_TABLE[5, row, 1] + + HELIO_LONG_TABLE[5, row, 2] * jme) + ) + + l_rad = (l0 + l1 * jme + l2 * jme**2 + l3 * jme**3 + l4 * jme**4 + + l5 * jme**5)/10**8 + l = np.rad2deg(l_rad) + return l % 360 + + +@jcompile('float64(float64)', nopython=True) +def heliocentric_latitude(jme): + b0 = 0.0 + b1 = 0.0 + for row in range(HELIO_LAT_TABLE.shape[1]): + b0 += (HELIO_LAT_TABLE[0, row, 0] + * np.cos(HELIO_LAT_TABLE[0, row, 1] + + HELIO_LAT_TABLE[0, row, 2] * jme) + ) + b1 += (HELIO_LAT_TABLE[1, row, 0] + * np.cos(HELIO_LAT_TABLE[1, row, 1] + + HELIO_LAT_TABLE[1, row, 2] * jme) + ) + + b_rad = (b0 + b1 * jme)/10**8 + b = np.rad2deg(b_rad) + return b + + +@jcompile('float64(float64)', nopython=True) +def heliocentric_radius_vector(jme): + r0 = 0.0 + r1 = 0.0 + r2 = 0.0 + r3 = 0.0 + r4 = 0.0 + for row in range(HELIO_RADIUS_TABLE.shape[1]): + r0 += (HELIO_RADIUS_TABLE[0, row, 0] + * np.cos(HELIO_RADIUS_TABLE[0, row, 1] + + HELIO_RADIUS_TABLE[0, row, 2] * jme) + ) + r1 += (HELIO_RADIUS_TABLE[1, row, 0] + * np.cos(HELIO_RADIUS_TABLE[1, row, 1] + + HELIO_RADIUS_TABLE[1, row, 2] * jme) + ) + r2 += (HELIO_RADIUS_TABLE[2, row, 0] + * np.cos(HELIO_RADIUS_TABLE[2, row, 1] + + HELIO_RADIUS_TABLE[2, row, 2] * jme) + ) + r3 += (HELIO_RADIUS_TABLE[3, row, 0] + * np.cos(HELIO_RADIUS_TABLE[3, row, 1] + + HELIO_RADIUS_TABLE[3, row, 2] * jme) + ) + r4 += (HELIO_RADIUS_TABLE[4, row, 0] + * np.cos(HELIO_RADIUS_TABLE[4, row, 1] + + HELIO_RADIUS_TABLE[4, row, 2] * jme) + ) + + r = (r0 + r1 * jme + r2 * jme**2 + r3 * jme**3 + r4 * jme**4)/10**8 + return r + + +@jcompile('float64(float64)', nopython=True) +def geocentric_longitude(heliocentric_longitude): + theta = heliocentric_longitude + 180.0 + return theta % 360 + + +@jcompile('float64(float64)', nopython=True) +def geocentric_latitude(heliocentric_latitude): + beta = -1.0*heliocentric_latitude + return beta + + +@jcompile('float64(float64)', nopython=True) +def mean_elongation(julian_ephemeris_century): + x0 = (297.85036 + + 445267.111480 * julian_ephemeris_century + - 0.0019142 * julian_ephemeris_century**2 + + julian_ephemeris_century**3 / 189474) + return x0 + + +@jcompile('float64(float64)', nopython=True) +def mean_anomaly_sun(julian_ephemeris_century): + x1 = (357.52772 + + 35999.050340 * julian_ephemeris_century + - 0.0001603 * julian_ephemeris_century**2 + - julian_ephemeris_century**3 / 300000) + return x1 + + +@jcompile('float64(float64)', nopython=True) +def mean_anomaly_moon(julian_ephemeris_century): + x2 = (134.96298 + + 477198.867398 * julian_ephemeris_century + + 0.0086972 * julian_ephemeris_century**2 + + julian_ephemeris_century**3 / 56250) + return x2 + + +@jcompile('float64(float64)', nopython=True) +def moon_argument_latitude(julian_ephemeris_century): + x3 = (93.27191 + + 483202.017538 * julian_ephemeris_century + - 0.0036825 * julian_ephemeris_century**2 + + julian_ephemeris_century**3 / 327270) + return x3 + + +@jcompile('float64(float64)', nopython=True) +def moon_ascending_longitude(julian_ephemeris_century): + x4 = (125.04452 + - 1934.136261 * julian_ephemeris_century + + 0.0020708 * julian_ephemeris_century**2 + + julian_ephemeris_century**3 / 450000) + return x4 + + +@jcompile('float64(float64, float64, float64, float64, float64, float64)', + nopython=True) +def longitude_nutation(julian_ephemeris_century, x0, x1, x2, x3, x4): + delta_psi_sum = 0 + for row in range(NUTATION_YTERM_ARRAY.shape[0]): + a = NUTATION_ABCD_ARRAY[row,0] + b = NUTATION_ABCD_ARRAY[row,1] + argsin = (NUTATION_YTERM_ARRAY[row,0]*x0 + + NUTATION_YTERM_ARRAY[row,1]*x1 + + NUTATION_YTERM_ARRAY[row,2]*x2 + + NUTATION_YTERM_ARRAY[row,3]*x3 + + NUTATION_YTERM_ARRAY[row,4]*x4) + term = (a + b * julian_ephemeris_century) * np.sin(np.radians(argsin)) + delta_psi_sum += term + delta_psi = delta_psi_sum*1.0/36000000 + return delta_psi + + +@jcompile('float64(float64, float64, float64, float64, float64, float64)', + nopython=True) +def obliquity_nutation(julian_ephemeris_century, x0, x1, x2, x3, x4): + delta_eps_sum = 0.0 + for row in range(NUTATION_YTERM_ARRAY.shape[0]): + c = NUTATION_ABCD_ARRAY[row,2] + d = NUTATION_ABCD_ARRAY[row,3] + argcos = (NUTATION_YTERM_ARRAY[row,0]*x0 + + NUTATION_YTERM_ARRAY[row,1]*x1 + + NUTATION_YTERM_ARRAY[row,2]*x2 + + NUTATION_YTERM_ARRAY[row,3]*x3 + + NUTATION_YTERM_ARRAY[row,4]*x4) + term = (c + d * julian_ephemeris_century) * np.cos(np.radians(argcos)) + delta_eps_sum += term + delta_eps = delta_eps_sum*1.0/36000000 + return delta_eps + + +@jcompile('float64(float64)', nopython=True) +def mean_ecliptic_obliquity(julian_ephemeris_millennium): + U = 1.0*julian_ephemeris_millennium/10 + e0 = (84381.448 - 4680.93 * U - 1.55 * U**2 + + 1999.25 * U**3 - 51.38 * U**4 - 249.67 * U**5 + - 39.05 * U**6 + 7.12 * U**7 + 27.87 * U**8 + + 5.79 * U**9 + 2.45 * U**10) + return e0 + + +@jcompile('float64(float64, float64)', nopython=True) +def true_ecliptic_obliquity(mean_ecliptic_obliquity, obliquity_nutation): + e0 = mean_ecliptic_obliquity + deleps = obliquity_nutation + e = e0*1.0/3600 + deleps + return e + + +@jcompile('float64(float64)', nopython=True) +def aberration_correction(earth_radius_vector): + deltau = -20.4898 / (3600 * earth_radius_vector) + return deltau + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def apparent_sun_longitude(geocentric_longitude, longitude_nutation, + aberration_correction): + lamd = geocentric_longitude + longitude_nutation + aberration_correction + return lamd + + +@jcompile('float64(float64, float64)', nopython=True) +def mean_sidereal_time(julian_day, julian_century): + v0 = (280.46061837 + 360.98564736629 * (julian_day - 2451545) + + 0.000387933 * julian_century**2 - julian_century**3 / 38710000) + return v0 % 360.0 + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def apparent_sidereal_time(mean_sidereal_time, longitude_nutation, + true_ecliptic_obliquity): + v= mean_sidereal_time + longitude_nutation * np.cos( + np.radians(true_ecliptic_obliquity)) + return v + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def geocentric_sun_right_ascension(apparent_sun_longitude, + true_ecliptic_obliquity, + geocentric_latitude): + num = (np.sin(np.radians(apparent_sun_longitude)) + * np.cos(np.radians(true_ecliptic_obliquity)) + - np.tan(np.radians(geocentric_latitude)) + * np.sin(np.radians(true_ecliptic_obliquity))) + alpha = np.degrees(np.arctan2(num, np.cos( + np.radians(apparent_sun_longitude)))) + return alpha % 360 + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def geocentric_sun_declination(apparent_sun_longitude, true_ecliptic_obliquity, + geocentric_latitude): + delta = np.degrees(np.arcsin(np.sin(np.radians(geocentric_latitude)) * + np.cos(np.radians(true_ecliptic_obliquity)) + + np.cos(np.radians(geocentric_latitude)) * + np.sin(np.radians(true_ecliptic_obliquity))* + np.sin(np.radians(apparent_sun_longitude)))) + return delta + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def local_hour_angle(apparent_sidereal_time, observer_longitude, + sun_right_ascension): + """Measured westward from south""" + H = apparent_sidereal_time + observer_longitude - sun_right_ascension + return H % 360 + + +@jcompile('float64(float64)', nopython=True) +def equatorial_horizontal_parallax(earth_radius_vector): + xi = 8.794 / (3600 * earth_radius_vector) + return xi + + +@jcompile('float64(float64)', nopython=True) +def uterm(observer_latitude): + u = np.arctan(0.99664719 * np.tan(np.radians(observer_latitude))) + return u + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def xterm(u, observer_latitude, observer_elevation): + x = (np.cos(u) + observer_elevation / 6378140 + * np.cos(np.radians(observer_latitude))) + return x + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def yterm(u, observer_latitude, observer_elevation): + y = (0.99664719 * np.sin(u) + observer_elevation / 6378140 + * np.sin(np.radians(observer_latitude))) + return y + + +@jcompile('float64(float64, float64,float64, float64)', nopython=True) +def parallax_sun_right_ascension(xterm, equatorial_horizontal_parallax, + local_hour_angle, geocentric_sun_declination): + num = (-xterm * np.sin(np.radians(equatorial_horizontal_parallax)) + * np.sin(np.radians(local_hour_angle))) + denom = (np.cos(np.radians(geocentric_sun_declination)) + - xterm * np.sin(np.radians(equatorial_horizontal_parallax)) + * np.cos(np.radians(local_hour_angle))) + delta_alpha = np.degrees(np.arctan2(num, denom)) + return delta_alpha + + +@jcompile('float64(float64, float64)', nopython=True) +def topocentric_sun_right_ascension(geocentric_sun_right_ascension, + parallax_sun_right_ascension): + alpha_prime = geocentric_sun_right_ascension + parallax_sun_right_ascension + return alpha_prime + + +@jcompile('float64(float64, float64, float64, float64, float64, float64)', + nopython=True) +def topocentric_sun_declination(geocentric_sun_declination, xterm, yterm, + equatorial_horizontal_parallax, + parallax_sun_right_ascension, + local_hour_angle): + num = ((np.sin(np.radians(geocentric_sun_declination)) - yterm + * np.sin(np.radians(equatorial_horizontal_parallax))) + * np.cos(np.radians(parallax_sun_right_ascension))) + denom = (np.cos(np.radians(geocentric_sun_declination)) - xterm + * np.sin(np.radians(equatorial_horizontal_parallax)) + * np.cos(np.radians(local_hour_angle))) + delta = np.degrees(np.arctan2(num, denom)) + return delta + + +@jcompile('float64(float64, float64)', nopython=True) +def topocentric_local_hour_angle(local_hour_angle, + parallax_sun_right_ascension): + H_prime = local_hour_angle - parallax_sun_right_ascension + return H_prime + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def topocentric_elevation_angle_without_atmosphere(observer_latitude, + topocentric_sun_declination, + topocentric_local_hour_angle): + e0 = np.degrees(np.arcsin( + np.sin(np.radians(observer_latitude)) + * np.sin(np.radians(topocentric_sun_declination)) + + np.cos(np.radians(observer_latitude)) + * np.cos(np.radians(topocentric_sun_declination)) + * np.cos(np.radians(topocentric_local_hour_angle)))) + return e0 + + +@jcompile('float64(float64, float64, float64, float64)', nopython=True) +def atmospheric_refraction_correction(local_pressure, local_temp, + topocentric_elevation_angle_wo_atmosphere, + atmos_refract + ): + switch = topocentric_elevation_angle_wo_atmosphere >= -1*(0.26667 + atmos_refract) + delta_e = ((local_pressure / 1010.0) * (283.0 / (273 + local_temp)) + * 1.02 / (60 * np.tan(np.radians( + topocentric_elevation_angle_wo_atmosphere + + 10.3 / (topocentric_elevation_angle_wo_atmosphere + + 5.11))))) * switch + return delta_e + + +@jcompile('float64(float64, float64)', nopython=True) +def topocentric_elevation_angle(topocentric_elevation_angle_without_atmosphere, + atmospheric_refraction_correction): + e = (topocentric_elevation_angle_without_atmosphere + + atmospheric_refraction_correction) + return e + + +@jcompile('float64(float64)', nopython=True) +def topocentric_zenith_angle(topocentric_elevation_angle): + theta = 90 - topocentric_elevation_angle + return theta + + +@jcompile('float64(float64, float64, float64)', nopython=True) +def topocentric_astronomers_azimuth(topocentric_local_hour_angle, + topocentric_sun_declination, + observer_latitude): + num = np.sin(np.radians(topocentric_local_hour_angle)) + denom = (np.cos(np.radians(topocentric_local_hour_angle)) + * np.sin(np.radians(observer_latitude)) + - np.tan(np.radians(topocentric_sun_declination)) + * np.cos(np.radians(observer_latitude))) + gamma = np.degrees(np.arctan2(num, denom)) + return gamma % 360 + + +@jcompile('float64(float64)', nopython=True) +def topocentric_azimuth_angle(topocentric_astronomers_azimuth): + phi = topocentric_astronomers_azimuth + 180 + return phi % 360 + + +@jcompile('void(float64[:], float64[:], float64[:,:])', nopython=True, + nogil=True) +def solar_position_loop(unixtime, loc_args, out): + """Loop through the time array and calculate the solar position""" + lat = loc_args[0] + lon = loc_args[1] + elev = loc_args[2] + pressure = loc_args[3] + temp = loc_args[4] + delta_t = loc_args[5] + atmos_refract = loc_args[6] + + for i in range(unixtime.shape[0]): + utime = unixtime[i] + jd = julian_day(utime) + jde = julian_ephemeris_day(jd, delta_t) + jc = julian_century(jd) + jce = julian_ephemeris_century(jde) + jme = julian_ephemeris_millennium(jce) + L = heliocentric_longitude(jme) + B = heliocentric_latitude(jme) + R = heliocentric_radius_vector(jme) + Theta = geocentric_longitude(L) + beta = geocentric_latitude(B) + x0 = mean_elongation(jce) + x1 = mean_anomaly_sun(jce) + x2 = mean_anomaly_moon(jce) + x3 = moon_argument_latitude(jce) + x4 = moon_ascending_longitude(jce) + delta_psi = longitude_nutation(jce, x0, x1, x2, x3, x4) + delta_epsilon = obliquity_nutation(jce, x0, x1, x2, x3, x4) + epsilon0 = mean_ecliptic_obliquity(jme) + epsilon = true_ecliptic_obliquity(epsilon0, delta_epsilon) + delta_tau = aberration_correction(R) + lamd = apparent_sun_longitude(Theta, delta_psi, delta_tau) + v0 = mean_sidereal_time(jd, jc) + v = apparent_sidereal_time(v0, delta_psi, epsilon) + alpha = geocentric_sun_right_ascension(lamd, epsilon, beta) + delta = geocentric_sun_declination(lamd, epsilon, beta) + H = local_hour_angle(v, lon, alpha) + xi = equatorial_horizontal_parallax(R) + u = uterm(lat) + x = xterm(u, lat, elev) + y = yterm(u, lat, elev) + delta_alpha = parallax_sun_right_ascension(x, xi, H, delta) + alpha_prime = topocentric_sun_right_ascension(alpha, delta_alpha) + delta_prime = topocentric_sun_declination(delta, x, y, xi, delta_alpha, + H) + H_prime = topocentric_local_hour_angle(H, delta_alpha) + e0 = topocentric_elevation_angle_without_atmosphere(lat, delta_prime, + H_prime) + delta_e = atmospheric_refraction_correction(pressure, temp, e0, + atmos_refract) + e = topocentric_elevation_angle(e0, delta_e) + theta = topocentric_zenith_angle(e) + theta0 = topocentric_zenith_angle(e0) + gamma = topocentric_astronomers_azimuth(H_prime, delta_prime, lat) + phi = topocentric_azimuth_angle(gamma) + out[0,i] = theta + out[1,i] = theta0 + out[2,i] = e + out[3,i] = e0 + out[4,i] = phi + + +def solar_position_numba(unixtime, lat, lon, elev, pressure, temp, delta_t, + atmos_refract, numthreads): + """Calculate the solar position using the numba compiled functions + and multiple threads. Very slow if functions are not numba compiled. + """ + loc_args = np.array([lat, lon, elev, pressure, temp, delta_t, + atmos_refract]) + ulength = unixtime.shape[0] + result = np.empty((5,ulength), dtype=np.float64) + + if ulength < numthreads: + pvl_logger.warning('The number of threads is more than the length of' + + ' the time array. Only using {} threads.'.format( + ulength)) + numthreads = ulength + + if numthreads <= 1: + pvl_logger.debug('Only using one thread for calculation') + solar_position_loop(unixtime, loc_args, result) + return result + + split0 = np.array_split(unixtime, numthreads) + split2 = np.array_split(result, numthreads, axis=1) + chunks = [ [a0, loc_args, split2[i]] for i, a0 in enumerate(split0)] + # Spawn one thread per chunk + threads = [threading.Thread(target=solar_position_loop, args=chunk) + for chunk in chunks] + for thread in threads: + thread.start() + for thread in threads: + thread.join() + return result + + +def solar_position_numpy(unixtime, lat, lon, elev, pressure, temp, delta_t, + atmos_refract, numthreads): + """Calculate the solar position assuming unixtime is a numpy array. Note + this function will not work if the solar position functions were + compiled with numba. + """ + + jd = julian_day(unixtime) + jde = julian_ephemeris_day(jd, delta_t) + jc = julian_century(jd) + jce = julian_ephemeris_century(jde) + jme = julian_ephemeris_millennium(jce) + L = heliocentric_longitude(jme) + B = heliocentric_latitude(jme) + R = heliocentric_radius_vector(jme) + Theta = geocentric_longitude(L) + beta = geocentric_latitude(B) + x0 = mean_elongation(jce) + x1 = mean_anomaly_sun(jce) + x2 = mean_anomaly_moon(jce) + x3 = moon_argument_latitude(jce) + x4 = moon_ascending_longitude(jce) + delta_psi = longitude_nutation(jce, x0, x1, x2, x3, x4) + delta_epsilon = obliquity_nutation(jce, x0, x1, x2, x3, x4) + epsilon0 = mean_ecliptic_obliquity(jme) + epsilon = true_ecliptic_obliquity(epsilon0, delta_epsilon) + delta_tau = aberration_correction(R) + lamd = apparent_sun_longitude(Theta, delta_psi, delta_tau) + v0 = mean_sidereal_time(jd, jc) + v = apparent_sidereal_time(v0, delta_psi, epsilon) + alpha = geocentric_sun_right_ascension(lamd, epsilon, beta) + delta = geocentric_sun_declination(lamd, epsilon, beta) + H = local_hour_angle(v, lon, alpha) + xi = equatorial_horizontal_parallax(R) + u = uterm(lat) + x = xterm(u, lat, elev) + y = yterm(u, lat, elev) + delta_alpha = parallax_sun_right_ascension(x, xi, H, delta) + alpha_prime = topocentric_sun_right_ascension(alpha, delta_alpha) + delta_prime = topocentric_sun_declination(delta, x, y, xi, delta_alpha, H) + H_prime = topocentric_local_hour_angle(H, delta_alpha) + e0 = topocentric_elevation_angle_without_atmosphere(lat, delta_prime, + H_prime) + delta_e = atmospheric_refraction_correction(pressure, temp, e0, atmos_refract) + e = topocentric_elevation_angle(e0, delta_e) + theta = topocentric_zenith_angle(e) + theta0 = topocentric_zenith_angle(e0) + gamma = topocentric_astronomers_azimuth(H_prime, delta_prime, lat) + phi = topocentric_azimuth_angle(gamma) + return theta, theta0, e, e0, phi + + +def solar_position(unixtime, lat, lon, elev, pressure, temp, delta_t, + atmos_refract, numthreads=8,): + + """ + Calculate the solar position using a translation of the + NREL SPA code. + + If numba is installed, the functions can be compiled + and the code runs quickly. If not, the functions + still evaluate but use numpy instead + + Parameters + ---------- + unixtime : numpy array + Array of unix/epoch timestamps to calculate solar position for. + Unixtime is the number of seconds since Jan. 1, 1970 00:00:00 UTC. + A pandas.DatetimeIndex is easily converted using .astype(int)/10**9 + lat : float + Latitude to calculate solar position for + lon : float + Longitude to calculate solar position for + elev : float + Elevation of location in meters + pressure : int or float + avg. yearly pressure at location in Pascals; + used for atmospheric correction + temp : int or float + avg. yearly temperature at location in + degrees C; used for atmospheric correction + numthreads: int, optional + Number of threads to use for computation if numba>=0.17 + is installed. + delta_t : float, optional + Difference between terrestrial time and UT1. + By default, use USNO historical data and predictions + + Returns + ------- + Numpy Array with elements: + apparent zenith, + zenith, + elevation, + apparent_elevation, + azimuth + + 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. + """ + if USE_NUMBA: + do_calc = solar_position_numba + else: + do_calc = solar_position_numpy + + + result = do_calc(unixtime, lat, lon, elev, pressure, + temp, delta_t, atmos_refract, numthreads) + + if not isinstance(result, np.ndarray): + try: + result = np.array(result) + except: + pass + + return result diff --git a/pvlib/spa_c_files/cspa_py.pxd b/pvlib/spa_c_files/cspa_py.pxd index 7ee51b91f7..bcf533493f 100644 --- a/pvlib/spa_c_files/cspa_py.pxd +++ b/pvlib/spa_c_files/cspa_py.pxd @@ -26,6 +26,8 @@ cdef extern from "spa.h": int function + double e0 + double e double zenith double azimuth_astro double azimuth diff --git a/pvlib/spa_c_files/spa_py.c b/pvlib/spa_c_files/spa_py.c index 12d21e214a..28f9127762 100644 --- a/pvlib/spa_c_files/spa_py.c +++ b/pvlib/spa_c_files/spa_py.c @@ -1,4 +1,4 @@ -/* Generated by Cython 0.21 */ +/* Generated by Cython 0.21.2 */ #define PY_SSIZE_T_CLEAN #ifndef CYTHON_USE_PYLONG_INTERNALS @@ -19,7 +19,7 @@ #elif PY_VERSION_HEX < 0x02060000 || (0x03000000 <= PY_VERSION_HEX && PY_VERSION_HEX < 0x03020000) #error Cython requires Python 2.6+ or Python 3.2+. #else -#define CYTHON_ABI "0_21" +#define CYTHON_ABI "0_21_2" #include #ifndef offsetof #define offsetof(type, member) ( (size_t) & ((type*)0) -> member ) @@ -73,8 +73,6 @@ #if PY_MAJOR_VERSION >= 3 #define Py_TPFLAGS_CHECKTYPES 0 #define Py_TPFLAGS_HAVE_INDEX 0 -#endif -#if PY_MAJOR_VERSION >= 3 #define Py_TPFLAGS_HAVE_NEWBUFFER 0 #endif #if PY_VERSION_HEX < 0x030400a1 && !defined(Py_TPFLAGS_HAVE_FINALIZE) @@ -101,10 +99,12 @@ #if CYTHON_COMPILING_IN_PYPY #define __Pyx_PyUnicode_Concat(a, b) PyNumber_Add(a, b) #define __Pyx_PyUnicode_ConcatSafe(a, b) PyNumber_Add(a, b) + #define __Pyx_PyFrozenSet_Size(s) PyObject_Size(s) #else #define __Pyx_PyUnicode_Concat(a, b) PyUnicode_Concat(a, b) #define __Pyx_PyUnicode_ConcatSafe(a, b) ((unlikely((a) == Py_None) || unlikely((b) == Py_None)) ? \ PyNumber_Add(a, b) : __Pyx_PyUnicode_Concat(a, b)) + #define __Pyx_PyFrozenSet_Size(s) PySet_Size(s) #endif #define __Pyx_PyString_FormatSafe(a, b) ((unlikely((a) == Py_None)) ? PyNumber_Remainder(a, b) : __Pyx_PyString_Format(a, b)) #define __Pyx_PyUnicode_FormatSafe(a, b) ((unlikely((a) == Py_None)) ? PyNumber_Remainder(a, b) : PyUnicode_Format(a, b)) @@ -151,6 +151,11 @@ #if PY_MAJOR_VERSION >= 3 #define PyBoolObject PyLongObject #endif +#if PY_MAJOR_VERSION >= 3 && CYTHON_COMPILING_IN_PYPY + #ifndef PyUnicode_InternFromString + #define PyUnicode_InternFromString(s) PyUnicode_FromString(s) + #endif +#endif #if PY_VERSION_HEX < 0x030200A4 typedef long Py_hash_t; #define __Pyx_PyInt_FromHash_t PyInt_FromLong @@ -160,7 +165,9 @@ #define __Pyx_PyInt_AsHash_t PyInt_AsSsize_t #endif #if PY_MAJOR_VERSION >= 3 - #define PyMethod_New(func, self, klass) ((self) ? PyMethod_New(func, self) : PyInstanceMethod_New(func)) + #define __Pyx_PyMethod_New(func, self, klass) ((self) ? PyMethod_New(func, self) : PyInstanceMethod_New(func)) +#else + #define __Pyx_PyMethod_New(func, self, klass) PyMethod_New(func, self, klass) #endif #ifndef CYTHON_INLINE #if defined(__GNUC__) @@ -522,7 +529,9 @@ static int __Pyx_InitStrings(__Pyx_StringTabEntry *t); int __pyx_module_is_main_pvlib__spa_c_files__spa_py = 0; /* Implementation of 'pvlib.spa_c_files.spa_py' */ -static PyObject *__pyx_pf_5pvlib_11spa_c_files_6spa_py_spa_calc(CYTHON_UNUSED PyObject *__pyx_self, PyObject *__pyx_v_year, PyObject *__pyx_v_month, PyObject *__pyx_v_day, PyObject *__pyx_v_hour, PyObject *__pyx_v_minute, PyObject *__pyx_v_second, PyObject *__pyx_v_timezone, PyObject *__pyx_v_latitude, PyObject *__pyx_v_longitude, PyObject *__pyx_v_elevation); /* proto */ +static PyObject *__pyx_pf_5pvlib_11spa_c_files_6spa_py_spa_calc(CYTHON_UNUSED PyObject *__pyx_self, PyObject *__pyx_v_year, PyObject *__pyx_v_month, PyObject *__pyx_v_day, PyObject *__pyx_v_hour, PyObject *__pyx_v_minute, PyObject *__pyx_v_second, PyObject *__pyx_v_timezone, PyObject *__pyx_v_latitude, PyObject *__pyx_v_longitude, PyObject *__pyx_v_elevation, PyObject *__pyx_v_pressure, PyObject *__pyx_v_temperature, PyObject *__pyx_v_delta_t); /* proto */ +static char __pyx_k_e[] = "e"; +static char __pyx_k_e0[] = "e0"; static char __pyx_k_day[] = "day"; static char __pyx_k_err[] = "err"; static char __pyx_k_spa[] = "spa"; @@ -554,7 +563,7 @@ static char __pyx_k_azm_rotation[] = "azm_rotation"; static char __pyx_k_atmos_refract[] = "atmos_refract"; static char __pyx_k_azimuth_astro[] = "azimuth_astro"; static char __pyx_k_pvlib_spa_c_files_spa_py[] = "pvlib.spa_c_files.spa_py"; -static char __pyx_k_home_will_git_repos_pvlib_pytho[] = "/home/will/git_repos/pvlib-python/pvlib/spa_c_files/spa_py.pyx"; +static char __pyx_k_home_tony_git_repos_pvlib_pytho[] = "/home/tony/git_repos/pvlib-python/pvlib/spa_c_files/spa_py.pyx"; static PyObject *__pyx_n_s_atmos_refract; static PyObject *__pyx_n_s_azimuth; static PyObject *__pyx_n_s_azimuth_astro; @@ -562,10 +571,12 @@ static PyObject *__pyx_n_s_azm_rotation; static PyObject *__pyx_n_s_day; static PyObject *__pyx_n_s_delta_t; static PyObject *__pyx_n_s_delta_ut1; +static PyObject *__pyx_n_s_e; +static PyObject *__pyx_n_s_e0; static PyObject *__pyx_n_s_elevation; static PyObject *__pyx_n_s_err; static PyObject *__pyx_n_s_function; -static PyObject *__pyx_kp_s_home_will_git_repos_pvlib_pytho; +static PyObject *__pyx_kp_s_home_tony_git_repos_pvlib_pytho; static PyObject *__pyx_n_s_hour; static PyObject *__pyx_n_s_incidence; static PyObject *__pyx_n_s_latitude; @@ -593,7 +604,7 @@ static PyObject *__pyx_codeobj__2; /* "pvlib/spa_c_files/spa_py.pyx":3 * cimport cspa_py * - * def spa_calc(year,month,day,hour,minute,second,timezone,latitude,longitude,elevation): # <<<<<<<<<<<<<< + * def spa_calc(year,month,day,hour,minute,second,timezone,latitude,longitude,elevation,pressure,temperature,delta_t): # <<<<<<<<<<<<<< * cdef cspa_py.spa_data spa * */ @@ -612,6 +623,9 @@ static PyObject *__pyx_pw_5pvlib_11spa_c_files_6spa_py_1spa_calc(PyObject *__pyx PyObject *__pyx_v_latitude = 0; PyObject *__pyx_v_longitude = 0; PyObject *__pyx_v_elevation = 0; + PyObject *__pyx_v_pressure = 0; + PyObject *__pyx_v_temperature = 0; + PyObject *__pyx_v_delta_t = 0; int __pyx_lineno = 0; const char *__pyx_filename = NULL; int __pyx_clineno = 0; @@ -619,12 +633,15 @@ static PyObject *__pyx_pw_5pvlib_11spa_c_files_6spa_py_1spa_calc(PyObject *__pyx __Pyx_RefNannyDeclarations __Pyx_RefNannySetupContext("spa_calc (wrapper)", 0); { - static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_year,&__pyx_n_s_month,&__pyx_n_s_day,&__pyx_n_s_hour,&__pyx_n_s_minute,&__pyx_n_s_second,&__pyx_n_s_timezone,&__pyx_n_s_latitude,&__pyx_n_s_longitude,&__pyx_n_s_elevation,0}; - PyObject* values[10] = {0,0,0,0,0,0,0,0,0,0}; + static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_year,&__pyx_n_s_month,&__pyx_n_s_day,&__pyx_n_s_hour,&__pyx_n_s_minute,&__pyx_n_s_second,&__pyx_n_s_timezone,&__pyx_n_s_latitude,&__pyx_n_s_longitude,&__pyx_n_s_elevation,&__pyx_n_s_pressure,&__pyx_n_s_temperature,&__pyx_n_s_delta_t,0}; + PyObject* values[13] = {0,0,0,0,0,0,0,0,0,0,0,0,0}; if (unlikely(__pyx_kwds)) { Py_ssize_t kw_args; const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args); switch (pos_args) { + case 13: values[12] = PyTuple_GET_ITEM(__pyx_args, 12); + case 12: values[11] = PyTuple_GET_ITEM(__pyx_args, 11); + case 11: values[10] = PyTuple_GET_ITEM(__pyx_args, 10); case 10: values[9] = PyTuple_GET_ITEM(__pyx_args, 9); case 9: values[8] = PyTuple_GET_ITEM(__pyx_args, 8); case 8: values[7] = PyTuple_GET_ITEM(__pyx_args, 7); @@ -646,53 +663,68 @@ static PyObject *__pyx_pw_5pvlib_11spa_c_files_6spa_py_1spa_calc(PyObject *__pyx case 1: if (likely((values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_month)) != 0)) kw_args--; else { - __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 10, 10, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 13, 13, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} } case 2: if (likely((values[2] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_day)) != 0)) kw_args--; else { - __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 10, 10, 2); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 13, 13, 2); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} } case 3: if (likely((values[3] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_hour)) != 0)) kw_args--; else { - __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 10, 10, 3); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 13, 13, 3); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} } case 4: if (likely((values[4] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_minute)) != 0)) kw_args--; else { - __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 10, 10, 4); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 13, 13, 4); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} } case 5: if (likely((values[5] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_second)) != 0)) kw_args--; else { - __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 10, 10, 5); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 13, 13, 5); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} } case 6: if (likely((values[6] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_timezone)) != 0)) kw_args--; else { - __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 10, 10, 6); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 13, 13, 6); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} } case 7: if (likely((values[7] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_latitude)) != 0)) kw_args--; else { - __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 10, 10, 7); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 13, 13, 7); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} } case 8: if (likely((values[8] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_longitude)) != 0)) kw_args--; else { - __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 10, 10, 8); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 13, 13, 8); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} } case 9: if (likely((values[9] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_elevation)) != 0)) kw_args--; else { - __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 10, 10, 9); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 13, 13, 9); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + } + case 10: + if (likely((values[10] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_pressure)) != 0)) kw_args--; + else { + __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 13, 13, 10); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + } + case 11: + if (likely((values[11] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_temperature)) != 0)) kw_args--; + else { + __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 13, 13, 11); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + } + case 12: + if (likely((values[12] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_delta_t)) != 0)) kw_args--; + else { + __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 13, 13, 12); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} } } if (unlikely(kw_args > 0)) { if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "spa_calc") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} } - } else if (PyTuple_GET_SIZE(__pyx_args) != 10) { + } else if (PyTuple_GET_SIZE(__pyx_args) != 13) { goto __pyx_L5_argtuple_error; } else { values[0] = PyTuple_GET_ITEM(__pyx_args, 0); @@ -705,6 +737,9 @@ static PyObject *__pyx_pw_5pvlib_11spa_c_files_6spa_py_1spa_calc(PyObject *__pyx values[7] = PyTuple_GET_ITEM(__pyx_args, 7); values[8] = PyTuple_GET_ITEM(__pyx_args, 8); values[9] = PyTuple_GET_ITEM(__pyx_args, 9); + values[10] = PyTuple_GET_ITEM(__pyx_args, 10); + values[11] = PyTuple_GET_ITEM(__pyx_args, 11); + values[12] = PyTuple_GET_ITEM(__pyx_args, 12); } __pyx_v_year = values[0]; __pyx_v_month = values[1]; @@ -716,23 +751,26 @@ static PyObject *__pyx_pw_5pvlib_11spa_c_files_6spa_py_1spa_calc(PyObject *__pyx __pyx_v_latitude = values[7]; __pyx_v_longitude = values[8]; __pyx_v_elevation = values[9]; + __pyx_v_pressure = values[10]; + __pyx_v_temperature = values[11]; + __pyx_v_delta_t = values[12]; } goto __pyx_L4_argument_unpacking_done; __pyx_L5_argtuple_error:; - __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 10, 10, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + __Pyx_RaiseArgtupleInvalid("spa_calc", 1, 13, 13, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L3_error;} __pyx_L3_error:; __Pyx_AddTraceback("pvlib.spa_c_files.spa_py.spa_calc", __pyx_clineno, __pyx_lineno, __pyx_filename); __Pyx_RefNannyFinishContext(); return NULL; __pyx_L4_argument_unpacking_done:; - __pyx_r = __pyx_pf_5pvlib_11spa_c_files_6spa_py_spa_calc(__pyx_self, __pyx_v_year, __pyx_v_month, __pyx_v_day, __pyx_v_hour, __pyx_v_minute, __pyx_v_second, __pyx_v_timezone, __pyx_v_latitude, __pyx_v_longitude, __pyx_v_elevation); + __pyx_r = __pyx_pf_5pvlib_11spa_c_files_6spa_py_spa_calc(__pyx_self, __pyx_v_year, __pyx_v_month, __pyx_v_day, __pyx_v_hour, __pyx_v_minute, __pyx_v_second, __pyx_v_timezone, __pyx_v_latitude, __pyx_v_longitude, __pyx_v_elevation, __pyx_v_pressure, __pyx_v_temperature, __pyx_v_delta_t); /* function exit code */ __Pyx_RefNannyFinishContext(); return __pyx_r; } -static PyObject *__pyx_pf_5pvlib_11spa_c_files_6spa_py_spa_calc(CYTHON_UNUSED PyObject *__pyx_self, PyObject *__pyx_v_year, PyObject *__pyx_v_month, PyObject *__pyx_v_day, PyObject *__pyx_v_hour, PyObject *__pyx_v_minute, PyObject *__pyx_v_second, PyObject *__pyx_v_timezone, PyObject *__pyx_v_latitude, PyObject *__pyx_v_longitude, PyObject *__pyx_v_elevation) { +static PyObject *__pyx_pf_5pvlib_11spa_c_files_6spa_py_spa_calc(CYTHON_UNUSED PyObject *__pyx_self, PyObject *__pyx_v_year, PyObject *__pyx_v_month, PyObject *__pyx_v_day, PyObject *__pyx_v_hour, PyObject *__pyx_v_minute, PyObject *__pyx_v_second, PyObject *__pyx_v_timezone, PyObject *__pyx_v_latitude, PyObject *__pyx_v_longitude, PyObject *__pyx_v_elevation, PyObject *__pyx_v_pressure, PyObject *__pyx_v_temperature, PyObject *__pyx_v_delta_t) { spa_data __pyx_v_spa; CYTHON_UNUSED int __pyx_v_err; PyObject *__pyx_r = NULL; @@ -810,7 +848,7 @@ static PyObject *__pyx_pf_5pvlib_11spa_c_files_6spa_py_spa_calc(CYTHON_UNUSED Py * spa.second = second * spa.timezone = timezone # <<<<<<<<<<<<<< * spa.delta_ut1 = 0 - * spa.delta_t = 67 + * spa.delta_t = delta_t */ __pyx_t_2 = __pyx_PyFloat_AsDouble(__pyx_v_timezone); if (unlikely((__pyx_t_2 == (double)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 12; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __pyx_v_spa.timezone = __pyx_t_2; @@ -819,7 +857,7 @@ static PyObject *__pyx_pf_5pvlib_11spa_c_files_6spa_py_spa_calc(CYTHON_UNUSED Py * spa.second = second * spa.timezone = timezone * spa.delta_ut1 = 0 # <<<<<<<<<<<<<< - * spa.delta_t = 67 + * spa.delta_t = delta_t * spa.longitude = longitude */ __pyx_v_spa.delta_ut1 = 0.0; @@ -827,15 +865,16 @@ static PyObject *__pyx_pf_5pvlib_11spa_c_files_6spa_py_spa_calc(CYTHON_UNUSED Py /* "pvlib/spa_c_files/spa_py.pyx":14 * spa.timezone = timezone * spa.delta_ut1 = 0 - * spa.delta_t = 67 # <<<<<<<<<<<<<< + * spa.delta_t = delta_t # <<<<<<<<<<<<<< * spa.longitude = longitude * spa.latitude = latitude */ - __pyx_v_spa.delta_t = 67.0; + __pyx_t_2 = __pyx_PyFloat_AsDouble(__pyx_v_delta_t); if (unlikely((__pyx_t_2 == (double)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 14; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_v_spa.delta_t = __pyx_t_2; /* "pvlib/spa_c_files/spa_py.pyx":15 * spa.delta_ut1 = 0 - * spa.delta_t = 67 + * spa.delta_t = delta_t * spa.longitude = longitude # <<<<<<<<<<<<<< * spa.latitude = latitude * spa.elevation = elevation @@ -844,11 +883,11 @@ static PyObject *__pyx_pf_5pvlib_11spa_c_files_6spa_py_spa_calc(CYTHON_UNUSED Py __pyx_v_spa.longitude = __pyx_t_2; /* "pvlib/spa_c_files/spa_py.pyx":16 - * spa.delta_t = 67 + * spa.delta_t = delta_t * spa.longitude = longitude * spa.latitude = latitude # <<<<<<<<<<<<<< * spa.elevation = elevation - * spa.pressure = 820 + * spa.pressure = pressure */ __pyx_t_2 = __pyx_PyFloat_AsDouble(__pyx_v_latitude); if (unlikely((__pyx_t_2 == (double)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 16; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __pyx_v_spa.latitude = __pyx_t_2; @@ -857,8 +896,8 @@ static PyObject *__pyx_pf_5pvlib_11spa_c_files_6spa_py_spa_calc(CYTHON_UNUSED Py * spa.longitude = longitude * spa.latitude = latitude * spa.elevation = elevation # <<<<<<<<<<<<<< - * spa.pressure = 820 - * spa.temperature = 11 + * spa.pressure = pressure + * spa.temperature = temperature */ __pyx_t_2 = __pyx_PyFloat_AsDouble(__pyx_v_elevation); if (unlikely((__pyx_t_2 == (double)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 17; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __pyx_v_spa.elevation = __pyx_t_2; @@ -866,24 +905,26 @@ static PyObject *__pyx_pf_5pvlib_11spa_c_files_6spa_py_spa_calc(CYTHON_UNUSED Py /* "pvlib/spa_c_files/spa_py.pyx":18 * spa.latitude = latitude * spa.elevation = elevation - * spa.pressure = 820 # <<<<<<<<<<<<<< - * spa.temperature = 11 + * spa.pressure = pressure # <<<<<<<<<<<<<< + * spa.temperature = temperature * spa.slope = 30 */ - __pyx_v_spa.pressure = 820.0; + __pyx_t_2 = __pyx_PyFloat_AsDouble(__pyx_v_pressure); if (unlikely((__pyx_t_2 == (double)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 18; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_v_spa.pressure = __pyx_t_2; /* "pvlib/spa_c_files/spa_py.pyx":19 * spa.elevation = elevation - * spa.pressure = 820 - * spa.temperature = 11 # <<<<<<<<<<<<<< + * spa.pressure = pressure + * spa.temperature = temperature # <<<<<<<<<<<<<< * spa.slope = 30 * spa.azm_rotation = -10 */ - __pyx_v_spa.temperature = 11.0; + __pyx_t_2 = __pyx_PyFloat_AsDouble(__pyx_v_temperature); if (unlikely((__pyx_t_2 == (double)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 19; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_v_spa.temperature = __pyx_t_2; /* "pvlib/spa_c_files/spa_py.pyx":20 - * spa.pressure = 820 - * spa.temperature = 11 + * spa.pressure = pressure + * spa.temperature = temperature * spa.slope = 30 # <<<<<<<<<<<<<< * spa.azm_rotation = -10 * spa.atmos_refract = 0.5667 @@ -891,7 +932,7 @@ static PyObject *__pyx_pf_5pvlib_11spa_c_files_6spa_py_spa_calc(CYTHON_UNUSED Py __pyx_v_spa.slope = 30.0; /* "pvlib/spa_c_files/spa_py.pyx":21 - * spa.temperature = 11 + * spa.temperature = temperature * spa.slope = 30 * spa.azm_rotation = -10 # <<<<<<<<<<<<<< * spa.atmos_refract = 0.5667 @@ -934,7 +975,7 @@ static PyObject *__pyx_pf_5pvlib_11spa_c_files_6spa_py_spa_calc(CYTHON_UNUSED Py /* "pvlib/spa_c_files/spa_py.pyx":3 * cimport cspa_py * - * def spa_calc(year,month,day,hour,minute,second,timezone,latitude,longitude,elevation): # <<<<<<<<<<<<<< + * def spa_calc(year,month,day,hour,minute,second,timezone,latitude,longitude,elevation,pressure,temperature,delta_t): # <<<<<<<<<<<<<< * cdef cspa_py.spa_data spa * */ @@ -980,10 +1021,12 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = { {&__pyx_n_s_day, __pyx_k_day, sizeof(__pyx_k_day), 0, 0, 1, 1}, {&__pyx_n_s_delta_t, __pyx_k_delta_t, sizeof(__pyx_k_delta_t), 0, 0, 1, 1}, {&__pyx_n_s_delta_ut1, __pyx_k_delta_ut1, sizeof(__pyx_k_delta_ut1), 0, 0, 1, 1}, + {&__pyx_n_s_e, __pyx_k_e, sizeof(__pyx_k_e), 0, 0, 1, 1}, + {&__pyx_n_s_e0, __pyx_k_e0, sizeof(__pyx_k_e0), 0, 0, 1, 1}, {&__pyx_n_s_elevation, __pyx_k_elevation, sizeof(__pyx_k_elevation), 0, 0, 1, 1}, {&__pyx_n_s_err, __pyx_k_err, sizeof(__pyx_k_err), 0, 0, 1, 1}, {&__pyx_n_s_function, __pyx_k_function, sizeof(__pyx_k_function), 0, 0, 1, 1}, - {&__pyx_kp_s_home_will_git_repos_pvlib_pytho, __pyx_k_home_will_git_repos_pvlib_pytho, sizeof(__pyx_k_home_will_git_repos_pvlib_pytho), 0, 0, 1, 0}, + {&__pyx_kp_s_home_tony_git_repos_pvlib_pytho, __pyx_k_home_tony_git_repos_pvlib_pytho, sizeof(__pyx_k_home_tony_git_repos_pvlib_pytho), 0, 0, 1, 0}, {&__pyx_n_s_hour, __pyx_k_hour, sizeof(__pyx_k_hour), 0, 0, 1, 1}, {&__pyx_n_s_incidence, __pyx_k_incidence, sizeof(__pyx_k_incidence), 0, 0, 1, 1}, {&__pyx_n_s_latitude, __pyx_k_latitude, sizeof(__pyx_k_latitude), 0, 0, 1, 1}, @@ -1018,14 +1061,14 @@ static int __Pyx_InitCachedConstants(void) { /* "pvlib/spa_c_files/spa_py.pyx":3 * cimport cspa_py * - * def spa_calc(year,month,day,hour,minute,second,timezone,latitude,longitude,elevation): # <<<<<<<<<<<<<< + * def spa_calc(year,month,day,hour,minute,second,timezone,latitude,longitude,elevation,pressure,temperature,delta_t): # <<<<<<<<<<<<<< * cdef cspa_py.spa_data spa * */ - __pyx_tuple_ = PyTuple_Pack(12, __pyx_n_s_year, __pyx_n_s_month, __pyx_n_s_day, __pyx_n_s_hour, __pyx_n_s_minute, __pyx_n_s_second, __pyx_n_s_timezone, __pyx_n_s_latitude, __pyx_n_s_longitude, __pyx_n_s_elevation, __pyx_n_s_spa, __pyx_n_s_err); if (unlikely(!__pyx_tuple_)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_tuple_ = PyTuple_Pack(15, __pyx_n_s_year, __pyx_n_s_month, __pyx_n_s_day, __pyx_n_s_hour, __pyx_n_s_minute, __pyx_n_s_second, __pyx_n_s_timezone, __pyx_n_s_latitude, __pyx_n_s_longitude, __pyx_n_s_elevation, __pyx_n_s_pressure, __pyx_n_s_temperature, __pyx_n_s_delta_t, __pyx_n_s_spa, __pyx_n_s_err); if (unlikely(!__pyx_tuple_)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_GOTREF(__pyx_tuple_); __Pyx_GIVEREF(__pyx_tuple_); - __pyx_codeobj__2 = (PyObject*)__Pyx_PyCode_New(10, 0, 12, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple_, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_home_will_git_repos_pvlib_pytho, __pyx_n_s_spa_calc, 3, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_codeobj__2 = (PyObject*)__Pyx_PyCode_New(13, 0, 15, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple_, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_home_tony_git_repos_pvlib_pytho, __pyx_n_s_spa_calc, 3, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 3; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_RefNannyFinishContext(); return 0; __pyx_L1_error:; @@ -1128,7 +1171,7 @@ PyMODINIT_FUNC PyInit_spa_py(void) /* "pvlib/spa_c_files/spa_py.pyx":3 * cimport cspa_py * - * def spa_calc(year,month,day,hour,minute,second,timezone,latitude,longitude,elevation): # <<<<<<<<<<<<<< + * def spa_calc(year,month,day,hour,minute,second,timezone,latitude,longitude,elevation,pressure,temperature,delta_t): # <<<<<<<<<<<<<< * cdef cspa_py.spa_data spa * */ @@ -1140,7 +1183,7 @@ PyMODINIT_FUNC PyInit_spa_py(void) /* "pvlib/spa_c_files/spa_py.pyx":1 * cimport cspa_py # <<<<<<<<<<<<<< * - * def spa_calc(year,month,day,hour,minute,second,timezone,latitude,longitude,elevation): + * def spa_calc(year,month,day,hour,minute,second,timezone,latitude,longitude,elevation,pressure,temperature,delta_t): */ __pyx_t_1 = PyDict_New(); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 1; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_GOTREF(__pyx_t_1); @@ -1155,7 +1198,6 @@ PyMODINIT_FUNC PyInit_spa_py(void) if (__pyx_m) { if (__pyx_d) { __Pyx_AddTraceback("init pvlib.spa_c_files.spa_py", __pyx_clineno, __pyx_lineno, __pyx_filename); - Py_DECREF(__pyx_d); __pyx_d = 0; } Py_DECREF(__pyx_m); __pyx_m = 0; } else if (!PyErr_Occurred()) { @@ -1685,6 +1727,12 @@ static PyObject* __pyx_convert__to_py_spa_data(spa_data s) { member = __Pyx_PyInt_From_int(s.function); if (member == NULL) goto bad; if (PyDict_SetItem(res, __pyx_n_s_function, member) < 0) goto bad; Py_DECREF(member); + member = PyFloat_FromDouble(s.e0); if (member == NULL) goto bad; + if (PyDict_SetItem(res, __pyx_n_s_e0, member) < 0) goto bad; + Py_DECREF(member); + member = PyFloat_FromDouble(s.e); if (member == NULL) goto bad; + if (PyDict_SetItem(res, __pyx_n_s_e, member) < 0) goto bad; + Py_DECREF(member); member = PyFloat_FromDouble(s.zenith); if (member == NULL) goto bad; if (PyDict_SetItem(res, __pyx_n_s_zenith, member) < 0) goto bad; Py_DECREF(member); diff --git a/pvlib/spa_c_files/spa_py.pyx b/pvlib/spa_c_files/spa_py.pyx index 36456ff4a6..2359045a63 100644 --- a/pvlib/spa_c_files/spa_py.pyx +++ b/pvlib/spa_c_files/spa_py.pyx @@ -1,6 +1,6 @@ cimport cspa_py -def spa_calc(year,month,day,hour,minute,second,timezone,latitude,longitude,elevation): +def spa_calc(year,month,day,hour,minute,second,timezone,latitude,longitude,elevation,pressure,temperature,delta_t): cdef cspa_py.spa_data spa spa.year = year @@ -11,12 +11,12 @@ def spa_calc(year,month,day,hour,minute,second,timezone,latitude,longitude,eleva spa.second = second spa.timezone = timezone spa.delta_ut1 = 0 - spa.delta_t = 67 + spa.delta_t = delta_t spa.longitude = longitude spa.latitude = latitude spa.elevation = elevation - spa.pressure = 820 - spa.temperature = 11 + spa.pressure = pressure + spa.temperature = temperature spa.slope = 30 spa.azm_rotation = -10 spa.atmos_refract = 0.5667 diff --git a/pvlib/test/test_solarposition.py b/pvlib/test/test_solarposition.py index 1e67bf0a71..a2f087c5c7 100644 --- a/pvlib/test/test_solarposition.py +++ b/pvlib/test/test_solarposition.py @@ -32,25 +32,27 @@ def test_spa_physical(): times = pd.date_range(datetime.datetime(2003,10,17,12,30,30), periods=1, freq='D') try: - ephem_data = solarposition.spa(times, golden_mst).ix[0] + ephem_data = solarposition.spa(times, golden_mst, pressure=82000, + temperature=11).ix[0] except ImportError: raise SkipTest - assert_almost_equals(50.111622, ephem_data['zenith'], 6) + assert_almost_equals(39.872046, ephem_data['elevation'], 6) + assert_almost_equals(50.111622, ephem_data['apparent_zenith'], 6) assert_almost_equals(194.340241, ephem_data['azimuth'], 6) - assert_almost_equals(39.888378, ephem_data['elevation'], 6) - + assert_almost_equals(39.888378, ephem_data['apparent_elevation'], 6) def test_spa_physical_dst(): times = pd.date_range(datetime.datetime(2003,10,17,13,30,30), periods=1, freq='D') try: - ephem_data = solarposition.spa(times, golden).ix[0] + ephem_data = solarposition.spa(times, golden, pressure=82000, + temperature=11).ix[0] except ImportError: raise SkipTest - assert_almost_equals(50.111622, ephem_data['zenith'], 6) + assert_almost_equals(39.872046, ephem_data['elevation'], 6) + assert_almost_equals(50.111622, ephem_data['apparent_zenith'], 6) assert_almost_equals(194.340241, ephem_data['azimuth'], 6) - assert_almost_equals(39.888378, ephem_data['elevation'], 6) - + assert_almost_equals(39.888378, ephem_data['apparent_elevation'], 6) def test_spa_localization(): @@ -60,6 +62,73 @@ def test_spa_localization(): raise SkipTest +def test_spa_python_numpy_physical(): + times = pd.date_range(datetime.datetime(2003,10,17,12,30,30), periods=1, freq='D') + ephem_data = solarposition.spa_python(times, golden_mst, pressure=82000, + temperature=11, delta_t=67, + atmos_refract=0.5667, + how='numpy').ix[0] + assert_almost_equals(39.872046, ephem_data['elevation'], 6) + assert_almost_equals(50.111622, ephem_data['apparent_zenith'], 6) + assert_almost_equals(194.340241, ephem_data['azimuth'], 6) + assert_almost_equals(39.888378, ephem_data['apparent_elevation'], 6) + + +def test_spa_python_numpy_physical_dst(): + times = pd.date_range(datetime.datetime(2003,10,17,13,30,30), periods=1, freq='D') + ephem_data = solarposition.spa_python(times, golden, pressure=82000, + temperature=11, delta_t=67, + atmos_refract=0.5667, + how='numpy').ix[0] + assert_almost_equals(50.111622, ephem_data['apparent_zenith'], 6) + assert_almost_equals(194.340241, ephem_data['azimuth'], 6) + assert_almost_equals(39.888378, ephem_data['apparent_elevation'], 6) + + +def test_spa_python_numba_physical(): + try: + import numba + except ImportError: + raise SkipTest + vers = numba.__version__.split('.') + if int(vers[0] + vers[1]) < 17: + raise SkipTest + + times = pd.date_range(datetime.datetime(2003,10,17,12,30,30), periods=1, freq='D') + ephem_data = solarposition.spa_python(times, golden_mst, pressure=82000, + temperature=11, delta_t=67, + atmos_refract=0.5667, + how='numba', numthreads=1).ix[0] + assert_almost_equals(39.872046, ephem_data['elevation'], 6) + assert_almost_equals(50.111622, ephem_data['apparent_zenith'], 6) + assert_almost_equals(194.340241, ephem_data['azimuth'], 6) + assert_almost_equals(39.888378, ephem_data['apparent_elevation'], 6) + + +def test_spa_python_numba_physical_dst(): + try: + import numba + except ImportError: + raise SkipTest + vers = numba.__version__.split('.') + if int(vers[0] + vers[1]) < 17: + raise SkipTest + + times = pd.date_range(datetime.datetime(2003,10,17,13,30,30), periods=1, freq='D') + ephem_data = solarposition.spa_python(times, golden, pressure=82000, + temperature=11, delta_t=67, + atmos_refract=0.5667, + how='numba', numthreads=1).ix[0] + assert_almost_equals(50.111622, ephem_data['apparent_zenith'], 6) + assert_almost_equals(194.340241, ephem_data['azimuth'], 6) + assert_almost_equals(39.888378, ephem_data['apparent_elevation'], 6) + + +def test_spa_python_localization(): + assert_frame_equal(solarposition.spa_python(times, tus), + solarposition.spa_python(times_localized, tus)) + + def test_pyephem_physical(): times = pd.date_range(datetime.datetime(2003,10,17,12,30,30), periods=1, freq='D') ephem_data = solarposition.pyephem(times, golden_mst, pressure=82000, temperature=11).ix[0] diff --git a/pvlib/test/test_spa.py b/pvlib/test/test_spa.py new file mode 100644 index 0000000000..457c315268 --- /dev/null +++ b/pvlib/test/test_spa.py @@ -0,0 +1,279 @@ +import os +import logging +pvl_logger = logging.getLogger('pvlib') +try: + from importlib import reload +except ImportError: + try: + from imp import reload + except ImportError: + pass + +import numpy as np +import numpy.testing as npt +import pandas as pd + +import unittest +from nose.tools import raises, assert_almost_equals +from nose.plugins.skip import SkipTest + +from pvlib.location import Location + + +try: + from numba import __version__ as numba_version + numba_version_int = int(numba_version.split('.')[0] + + numba_version.split('.')[1]) +except ImportError: + numba_version_int = 0 + + +times = pd.date_range('2003-10-17 12:30:30', periods=1, freq='D').tz_localize('MST') +unixtimes = times.tz_convert('UTC').astype(int)*1.0/10**9 +lat = 39.742476 +lon = -105.1786 +elev = 1830.14 +pressure = 820 +temp = 11 +delta_t = 67.0 +atmos_refract= 0.5667 + +JD = 2452930.312847 +JC = 0.0379277986858 +JDE = 2452930.313623 +JCE = 0.037927819916852 +JME = 0.003792781991685 +L = 24.0182616917 +B = -0.0001011219 +R = 0.9965422974 +Theta = 204.0182616917 +beta = 0.0001011219 +X0 = 17185.861179 +X1 = 1722.893218 +X2 = 18234.075703 +X3 = 18420.071012 +X4 = 51.686951 +dPsi = -0.00399840 +dEpsilon = 0.00166657 +epsilon0 = 84379.672625 +epsilon = 23.440465 +dTau = -0.005711 +lamd = 204.0085519281 +v0 = 318.515579 +v = 318.511910 +alpha = 202.227408 +delta = -9.31434 +H = 11.10590 +xi = 0.002451 +dAlpha = -0.000369 +alpha_prime = 202.22704 +delta_prime = -9.316179 +H_prime = 11.10627 +e0 = 39.872046 +de = 0.016332 +e = 39.888378 +theta = 50.11162 +theta0 = 90 - e0 +Gamma = 14.340241 +Phi = 194.340241 + + +class SpaBase(object): + """Test functions common to numpy and numba spa""" + def test_julian_day_dt(self): + dt = times.tz_convert('UTC')[0] + year = dt.year + month = dt.month + day = dt.day + hour = dt.hour + minute = dt.minute + second = dt.second + microsecond = dt.microsecond + assert_almost_equals(JD, + self.spa.julian_day_dt(year, month, day, hour, + minute, second, microsecond), 6) + + def test_julian_ephemeris_day(self): + assert_almost_equals(JDE, self.spa.julian_ephemeris_day(JD, delta_t), 5) + + def test_julian_century(self): + assert_almost_equals(JC, self.spa.julian_century(JD), 6) + + def test_julian_ephemeris_century(self): + assert_almost_equals(JCE, self.spa.julian_ephemeris_century(JDE), 10) + + def test_julian_ephemeris_millenium(self): + assert_almost_equals(JME, self.spa.julian_ephemeris_millennium(JCE), 10) + + def test_heliocentric_longitude(self): + assert_almost_equals(L, self.spa.heliocentric_longitude(JME), 6) + + def test_heliocentric_latitude(self): + assert_almost_equals(B, self.spa.heliocentric_latitude(JME), 6) + + def test_heliocentric_radius_vector(self): + assert_almost_equals(R, self.spa.heliocentric_radius_vector(JME), 6) + + def test_geocentric_longitude(self): + assert_almost_equals(Theta, self.spa.geocentric_longitude(L), 6) + + def test_geocentric_latitude(self): + assert_almost_equals(beta, self.spa.geocentric_latitude(B), 6) + + def test_mean_elongation(self): + assert_almost_equals(X0, self.spa.mean_elongation(JCE), 5) + + def test_mean_anomaly_sun(self): + assert_almost_equals(X1, self.spa.mean_anomaly_sun(JCE), 5) + + def test_mean_anomaly_moon(self): + assert_almost_equals(X2, self.spa.mean_anomaly_moon(JCE), 5) + + def test_moon_argument_latitude(self): + assert_almost_equals(X3, self.spa.moon_argument_latitude(JCE), 5) + + def test_moon_ascending_longitude(self): + assert_almost_equals(X4, self.spa.moon_ascending_longitude(JCE), 6) + + def test_longitude_nutation(self): + assert_almost_equals(dPsi, self.spa.longitude_nutation(JCE, X0, X1, X2, + X3, X4), 6) + + def test_obliquity_nutation(self): + assert_almost_equals(dEpsilon, self.spa.obliquity_nutation(JCE, X0, X1, + X2, X3, X4), + 6) + + def test_mean_ecliptic_obliquity(self): + assert_almost_equals(epsilon0, self.spa.mean_ecliptic_obliquity(JME), 6) + + def test_true_ecliptic_obliquity(self): + assert_almost_equals(epsilon, self.spa.true_ecliptic_obliquity( + epsilon0, dEpsilon), 6) + + def test_aberration_correction(self): + assert_almost_equals(dTau, self.spa.aberration_correction(R), 6) + + def test_apparent_sun_longitude(self): + assert_almost_equals(lamd, self.spa.apparent_sun_longitude( + Theta, dPsi, dTau), 6) + + def test_mean_sidereal_time(self): + assert_almost_equals(v0, self.spa.mean_sidereal_time(JD, JC), 3) + + def test_apparent_sidereal_time(self): + assert_almost_equals(v, self.spa.apparent_sidereal_time( + v0, dPsi, epsilon), 5) + + def test_geocentric_sun_right_ascension(self): + assert_almost_equals(alpha, self.spa.geocentric_sun_right_ascension( + lamd, epsilon, beta), 6) + + def test_geocentric_sun_declination(self): + assert_almost_equals(delta, self.spa.geocentric_sun_declination( + lamd, epsilon, beta), 6) + + def test_local_hour_angle(self): + assert_almost_equals(H, self.spa.local_hour_angle(v, lon, alpha), 4) + + def test_equatorial_horizontal_parallax(self): + assert_almost_equals(xi, self.spa.equatorial_horizontal_parallax(R), 6) + + def test_parallax_sun_right_ascension(self): + u = self.spa.uterm(lat) + x = self.spa.xterm(u, lat, elev) + y = self.spa.yterm(u, lat, elev) + assert_almost_equals(dAlpha, self.spa.parallax_sun_right_ascension( + x, xi, H, delta), 4) + + def test_topocentric_sun_right_ascension(self): + assert_almost_equals(alpha_prime, + self.spa.topocentric_sun_right_ascension( + alpha, dAlpha), 5) + + def test_topocentric_sun_declination(self): + u = self.spa.uterm(lat) + x = self.spa.xterm(u, lat, elev) + y = self.spa.yterm(u, lat, elev) + assert_almost_equals(delta_prime, self.spa.topocentric_sun_declination( + delta, x, y, xi, dAlpha,H), 5) + + def test_topocentric_local_hour_angle(self): + assert_almost_equals(H_prime, self.spa.topocentric_local_hour_angle( + H, dAlpha), 5) + + def test_topocentric_elevation_angle_without_atmosphere(self): + assert_almost_equals( + e0, self.spa.topocentric_elevation_angle_without_atmosphere( + lat, delta_prime, H_prime), 6) + + def test_atmospheric_refraction_correction(self): + assert_almost_equals(de, self.spa.atmospheric_refraction_correction( + pressure, temp, e0, atmos_refract), 6) + + def test_topocentric_elevation_angle(self): + assert_almost_equals(e, self.spa.topocentric_elevation_angle(e0, de), 6) + + def test_topocentric_zenith_angle(self): + assert_almost_equals(theta, self.spa.topocentric_zenith_angle(e), 5) + + def test_topocentric_astronomers_azimuth(self): + assert_almost_equals(Gamma, self.spa.topocentric_astronomers_azimuth( + H_prime, delta_prime, lat), 5) + + def test_topocentric_azimuth_angle(self): + assert_almost_equals(Phi, self.spa.topocentric_azimuth_angle(Gamma), 5) + + def test_solar_position(self): + npt.assert_almost_equal( + np.array([[theta, theta0, e, e0, Phi]]).T, self.spa.solar_position( + unixtimes, lat, lon, elev, pressure, temp, delta_t, + atmos_refract), 5) + + +class NumpySpaTest(unittest.TestCase, SpaBase): + """Import spa without compiling to numba then run tests""" + @classmethod + def setUpClass(self): + os.environ['PVLIB_USE_NUMBA'] = '0' + import pvlib.spa as spa + spa = reload(spa) + self.spa = spa + + @classmethod + def tearDownClass(self): + del os.environ['PVLIB_USE_NUMBA'] + + def test_julian_day(self): + assert_almost_equals(JD, self.spa.julian_day(unixtimes)[0], 6) + + +@unittest.skipIf(numba_version_int < 17, + 'Numba not installed or version not >= 0.17.0') +class NumbaSpaTest(unittest.TestCase, SpaBase): + """Import spa, compiling to numba, and run tests""" + @classmethod + def setUpClass(self): + os.environ['PVLIB_USE_NUMBA'] = '1' + if numba_version_int >= 17: + import pvlib.spa as spa + spa = reload(spa) + self.spa = spa + + @classmethod + def tearDownClass(self): + del os.environ['PVLIB_USE_NUMBA'] + + def test_julian_day(self): + assert_almost_equals(JD, self.spa.julian_day(unixtimes[0]), 6) + + def test_solar_position_multithreaded(self): + result = np.array([theta, theta0, e, e0, Phi]) + nresult = np.array([result, result, result]).T + times = np.array([unixtimes[0], unixtimes[0], unixtimes[0]]) + npt.assert_almost_equal( + nresult + , self.spa.solar_position( + times, lat, lon, elev, pressure, temp, delta_t, + atmos_refract, numthreads=8), 5) + From 5cf27d2d2dedc3ad0c1a5c64603d079537b647f4 Mon Sep 17 00:00:00 2001 From: Tony Lorenzo Date: Wed, 8 Apr 2015 16:22:26 -0700 Subject: [PATCH 2/2] added sunrise/set/transit to python spa, removed pyephem dependency --- pvlib/solarposition.py | 230 ++++++++++++++++++++++++------- pvlib/spa.py | 208 ++++++++++++++++++++++++---- pvlib/test/test_solarposition.py | 47 ++++++- pvlib/test/test_spa.py | 103 +++++++++++++- setup.py | 1 - 5 files changed, 510 insertions(+), 79 deletions(-) diff --git a/pvlib/solarposition.py b/pvlib/solarposition.py index 24cf9c98e9..c73fadd51c 100644 --- a/pvlib/solarposition.py +++ b/pvlib/solarposition.py @@ -9,6 +9,7 @@ from __future__ import division import os +import warnings import logging pvl_logger = logging.getLogger('pvlib') import datetime as dt @@ -29,8 +30,8 @@ from pvlib.tools import localize_to_utc, datetime_to_djd, djd_to_datetime -def get_solarposition(time, location, method='pyephem', pressure=101325, - temperature=12): +def get_solarposition(time, location, method='nrel_numpy', pressure=101325, + temperature=12, **kwargs): """ A convenience wrapper for the solar position calculators. @@ -39,51 +40,63 @@ def get_solarposition(time, location, method='pyephem', pressure=101325, time : pandas.DatetimeIndex location : pvlib.Location object method : string - 'pyephem' uses the PyEphem package (default): :func:`pyephem` + 'pyephem' uses the PyEphem package: :func:`pyephem` - 'spa' or 'spa_c' uses the spa code: :func:`spa` + 'nrel_c' uses the NREL SPA C code [3]: :func:`spa_c` - 'spa_numpy' uses SPA code translated to python: :func: `spa_python` + 'nrel_numpy' uses an implementation of the NREL SPA algorithm + described in [1] (default): :func:`spa_python` - 'spa_numba' uses SPA code translated to python: :func: `spa_python` + 'nrel_numba' uses an implementation of the NREL SPA algorithm + described in [1], but also compiles the code first: :func:`spa_python` 'ephemeris' uses the pvlib ephemeris code: :func:`ephemeris` pressure : float Pascals. temperature : float Degrees C. + + Other keywords are passed to the underlying solar position function. + + 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: http://rredc.nrel.gov/solar/codesandalgorithms/spa/ """ method = method.lower() if isinstance(time, dt.datetime): time = pd.DatetimeIndex([time, ]) - if method == 'spa' or method == 'spa_c': - ephem_df = spa(time, location, pressure, temperature) - elif method == 'pyephem': - ephem_df = pyephem(time, location, pressure, temperature) - elif method == 'spa_numpy': + if method == 'nrel_c': + ephem_df = spa_c(time, location, pressure, temperature, **kwargs) + elif method == 'nrel_numba': ephem_df = spa_python(time, location, pressure, temperature, - how='numpy') - elif method == 'spa_numba': + how='numba', **kwargs) + elif method == 'nrel_numpy': ephem_df = spa_python(time, location, pressure, temperature, - how='numba') + how='numpy', **kwargs) + elif method == 'pyephem': + ephem_df = pyephem(time, location, pressure, temperature, **kwargs) elif method == 'ephemeris': - ephem_df = ephemeris(time, location, pressure, temperature) + ephem_df = ephemeris(time, location, pressure, temperature, **kwargs) else: raise ValueError('Invalid solar position method') return ephem_df -def spa(time, location, pressure=101325, temperature=12, delta_t=67.0, +def spa_c(time, location, pressure=101325, temperature=12, delta_t=67.0, raw_spa_output=False): ''' 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. + a README file which describes how the C code is wrapped in Python. + Due to license restrictions, the C code must be downloaded seperately + and used in accordance with it's license. Parameters ---------- @@ -113,6 +126,10 @@ def spa(time, location, pressure=101325, temperature=12, delta_t=67.0, ---------- NREL SPA code: http://rredc.nrel.gov/solar/codesandalgorithms/spa/ USNO delta T: http://www.usno.navy.mil/USNO/earth-orientation/eo-products/long-term + + See also + -------- + pyephem, spa_python, ephemeris ''' # Added by Rob Andrews (@Calama-Consulting), Calama Consulting, 2014 @@ -161,11 +178,42 @@ def spa(time, location, pressure=101325, temperature=12, delta_t=67.0, 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' + pvl_logger.debug('Reloading spa module without compiling') + spa = reload(spa) + 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' + pvl_logger.debug('Reloading spa module, compiling with numba') + spa = reload(spa) + del os.environ['PVLIB_USE_NUMBA'] + elif how != 'numba' and how != 'numpy': + raise ValueError("how must be either 'numba' or 'numpy'") + + return spa + + + def spa_python(time, location, pressure=101325, temperature=12, delta_t=None, atmos_refract=None, how='numpy', numthreads=4): """ - Calculate the solar position using a python translation of the - NREL SPA code. + Calculate the solar position using a python implementation of the + NREL SPA algorithm described in [1]. If numba is installed, the functions can be compiled to machine code and the function can be multithreaded. @@ -182,12 +230,14 @@ def spa_python(time, location, pressure=101325, temperature=12, delta_t=None, avg. yearly air temperature in degrees C. delta_t : float, optional Difference between terrestrial time and UT1. - By default, use USNO historical data and predictions + 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 - If numba >= 0.17.0 is installed, how='numba' will - compile the spa functions to machine code and - run them multithreaded. Otherwise, numpy calculations - are used. + 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 Number of threads to use if how == 'numba'. @@ -195,20 +245,27 @@ def spa_python(time, location, pressure=101325, temperature=12, delta_t=None, ------- DataFrame The DataFrame will have the following columns: - apparent_zenith, zenith, - apparent_elevation, elevation, - azimuth. + 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: http://www.usno.navy.mil/USNO/earth-orientation/eo-products/long-term + + See also + -------- + pyephem, spa_c, ephemeris """ # Added by Tony Lorenzo (@alorenzo175), University of Arizona, 2015 - from pvlib import spa pvl_logger.debug('Calculating solar position with spa_python code') lat = location.latitude @@ -226,26 +283,16 @@ def spa_python(time, location, pressure=101325, temperature=12, delta_t=None, unixtime = localize_to_utc(time, location).astype(int)/10**9 - using_numba = spa.USE_NUMBA - if how == 'numpy' and using_numba: - os.environ['PVLIB_USE_NUMBA'] = '0' - pvl_logger.debug('Reloading spa module without compiling') - spa = reload(spa) - del os.environ['PVLIB_USE_NUMBA'] - elif how == 'numba' and not using_numba: - os.environ['PVLIB_USE_NUMBA'] = '1' - pvl_logger.debug('Reloading spa module, compiling with numba') - spa = reload(spa) - del os.environ['PVLIB_USE_NUMBA'] - elif how != 'numba' and how != 'numpy': - raise ValueError("how must be either 'numba' or 'numpy'") - - app_zenith, zenith, app_elevation, elevation, azimuth = spa.solar_position( + spa = _spa_python_import(how) + + app_zenith, zenith, app_elevation, elevation, azimuth, eot = spa.solar_position( unixtime, lat, lon, elev, pressure, temperature, delta_t, atmos_refract, numthreads) + result = pd.DataFrame({'apparent_zenith':app_zenith, 'zenith':zenith, 'apparent_elevation':app_elevation, - 'elevation':elevation, 'azimuth':azimuth}, + 'elevation':elevation, 'azimuth':azimuth, + 'equation_of_time':eot}, index = time) try: @@ -254,6 +301,85 @@ def spa_python(time, location, pressure=101325, temperature=12, delta_t=None, result = result.tz_localize(location.tz) return result + + +def get_sun_rise_set_transit(time, location, how='numpy', delta_t=None, + numthreads=4): + """ + Calculate the sunrise, sunset, and sun transit times using the + NREL SPA algorithm 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 + Only the date part is used + location : pvlib.Location object + delta_t : float, optional + Difference between terrestrial time and UT1. + By default, use USNO historical data and predictions + how : str, optional + 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 + Number of threads to use if how == 'numba'. + + Returns + ------- + DataFrame + The DataFrame will have the following columns: + sunrise, sunset, 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 + + pvl_logger.debug('Calculating sunrise, set, transit with spa_python code') + + lat = location.latitude + lon = location.longitude + delta_t = delta_t or 67.0 + + if not isinstance(time, pd.DatetimeIndex): + try: + time = pd.DatetimeIndex(time) + except (TypeError, ValueError): + time = pd.DatetimeIndex([time,]) + + # must convert to midnight UTC on day of interest + utcday = pd.DatetimeIndex(time.date).tz_localize('UTC') + unixtime = utcday.astype(int)/10**9 + + spa = _spa_python_import(how) + + transit, sunrise, sunset = spa.transit_sunrise_sunset( + unixtime, lat, lon, delta_t, numthreads) + + # arrays are in seconds since epoch format, need to conver to timestamps + transit = pd.to_datetime(transit, unit='s', utc=True).tz_convert( + location.tz).tolist() + sunrise = pd.to_datetime(sunrise, unit='s', utc=True).tz_convert( + location.tz).tolist() + sunset = pd.to_datetime(sunset, unit='s', utc=True).tz_convert( + location.tz).tolist() + + result = pd.DataFrame({'transit':transit, + 'sunrise':sunrise, + 'sunset':sunset}, index=time) + + try: + result = result.tz_convert(location.tz) + except TypeError: + result = result.tz_localize(location.tz) + + return result def _ephem_setup(location, pressure, temperature): @@ -291,11 +417,18 @@ def pyephem(time, location, pressure=101325, temperature=12): apparent_elevation, elevation, apparent_azimuth, azimuth, apparent_zenith, zenith. + + See also + -------- + spa_python, spa_c, ephemeris + """ # Written by Will Holmgren (@wholmgren), University of Arizona, 2014 - - import ephem + try: + import ephem + except ImportError: + raise ImportError('PyEphem must be installed') pvl_logger.debug('using PyEphem to calculate solar position') @@ -382,7 +515,8 @@ def ephemeris(time, location, pressure=101325, temperature=12): See also -------- - pyephem, spa + pyephem, spa_c, spa_python + ''' # Added by Rob Andrews (@Calama-Consulting), Calama Consulting, 2014 @@ -594,3 +728,5 @@ def pyephem_earthsun_distance(time): earthsun.append(sun.earth_distance) return pd.Series(earthsun, index=time) + + diff --git a/pvlib/spa.py b/pvlib/spa.py index 0dac087588..7646075a3d 100644 --- a/pvlib/spa.py +++ b/pvlib/spa.py @@ -1,5 +1,5 @@ """ -Calculate the solar position using a translation of NREL SPA code either using +Calculate the solar position using the NREL SPA algorithm either using numpy arrays or compiling the code to machine language with numba. """ @@ -17,13 +17,20 @@ import numpy as np +# this block is a way to use an environment variable to switch between +# compiling the functions with numba or just use numpy def nocompile(*args, **kwargs): return lambda func: func if os.getenv('PVLIB_USE_NUMBA', '0') != '0': try: - from numba import vectorize, jit, float64, guvectorize, __version__ + from numba import jit, __version__ + except ImportError: + warnings.warn('Could not import numba, falling back to numpy calculation') + jcompile = nocompile + USE_NUMBA = False + else: major, minor = __version__.split('.')[:2] if int(major + minor) >= 17: # need at least numba >= 0.17.0 @@ -33,10 +40,6 @@ def nocompile(*args, **kwargs): warnings.warn('Numba version must be >= 0.17.0, falling back to numpy') jcompile = nocompile USE_NUMBA = False - except ImportError: - warnings.warn('Could not import numba, falling back to numpy calculation') - jcompile = nocompile - USE_NUMBA = False else: jcompile = nocompile USE_NUMBA = False @@ -442,7 +445,7 @@ def julian_day_dt(year, month, day, hour, minute, second, microsecond): @jcompile('float64(float64)', nopython=True) def julian_day(unixtime): - jd = unixtime*1.0/86400 + 2440587.5 + jd = unixtime * 1.0 / 86400 + 2440587.5 return jd @@ -818,9 +821,10 @@ def topocentric_elevation_angle_without_atmosphere(observer_latitude, @jcompile('float64(float64, float64, float64, float64)', nopython=True) def atmospheric_refraction_correction(local_pressure, local_temp, topocentric_elevation_angle_wo_atmosphere, - atmos_refract - ): - switch = topocentric_elevation_angle_wo_atmosphere >= -1*(0.26667 + atmos_refract) + atmos_refract): + # switch sets delta_e when the sun is below the horizon + switch = topocentric_elevation_angle_wo_atmosphere >= -1.0 * (0.26667 + + atmos_refract) delta_e = ((local_pressure / 1010.0) * (283.0 / (273 + local_temp)) * 1.02 / (60 * np.tan(np.radians( topocentric_elevation_angle_wo_atmosphere @@ -862,6 +866,32 @@ def topocentric_azimuth_angle(topocentric_astronomers_azimuth): return phi % 360 +@jcompile('float64(float64)', nopython=True) +def sun_mean_longitude(julian_ephemeris_millennium): + M = (280.4664567 + 360007.6982779 * julian_ephemeris_millennium + + 0.03032028 * julian_ephemeris_millennium**2 + + julian_ephemeris_millennium**3 / 49931 + - julian_ephemeris_millennium**4 / 15300 + - julian_ephemeris_millennium**5 / 2000000) + return M + + +@jcompile('float64(float64, float64, float64, float64)', nopython=True) +def equation_of_time(sun_mean_longitude, geocentric_sun_right_ascension, + longitude_nutation, true_ecliptic_obliquity): + E = (sun_mean_longitude - 0.0057183 - geocentric_sun_right_ascension + + longitude_nutation * np.cos(np.radians(true_ecliptic_obliquity))) + # limit between 0 and 360 + E = E % 360 + # convert to minutes + E *= 4 + greater = E > 20 + less = E < -20 + other = (E <= 20) & (E >= -20) + E = greater * (E - 1440) + less * (E + 1440) + other * E + return E + + @jcompile('void(float64[:], float64[:], float64[:,:])', nopython=True, nogil=True) def solar_position_loop(unixtime, loc_args, out): @@ -873,6 +903,7 @@ def solar_position_loop(unixtime, loc_args, out): temp = loc_args[4] delta_t = loc_args[5] atmos_refract = loc_args[6] + sst = loc_args[7] for i in range(unixtime.shape[0]): utime = unixtime[i] @@ -901,6 +932,13 @@ def solar_position_loop(unixtime, loc_args, out): v = apparent_sidereal_time(v0, delta_psi, epsilon) alpha = geocentric_sun_right_ascension(lamd, epsilon, beta) delta = geocentric_sun_declination(lamd, epsilon, beta) + if sst: + out[0,i] = v + out[1,i] = alpha + out[2,i] = delta + continue + m = sun_mean_longitude(jme) + eot = equation_of_time(m, alpha, delta_psi, epsilon) H = local_hour_angle(v, lon, alpha) xi = equatorial_horizontal_parallax(R) u = uterm(lat) @@ -925,17 +963,20 @@ def solar_position_loop(unixtime, loc_args, out): out[2,i] = e out[3,i] = e0 out[4,i] = phi - + out[5,i] = eot + def solar_position_numba(unixtime, lat, lon, elev, pressure, temp, delta_t, - atmos_refract, numthreads): + atmos_refract, numthreads, sst=False): """Calculate the solar position using the numba compiled functions and multiple threads. Very slow if functions are not numba compiled. """ loc_args = np.array([lat, lon, elev, pressure, temp, delta_t, - atmos_refract]) + atmos_refract, sst]) ulength = unixtime.shape[0] - result = np.empty((5,ulength), dtype=np.float64) + result = np.empty((6,ulength), dtype=np.float64) + if unixtime.dtype != np.float64: + unixtime = unixtime.astype(np.float64) if ulength < numthreads: pvl_logger.warning('The number of threads is more than the length of' + @@ -962,7 +1003,7 @@ def solar_position_numba(unixtime, lat, lon, elev, pressure, temp, delta_t, def solar_position_numpy(unixtime, lat, lon, elev, pressure, temp, delta_t, - atmos_refract, numthreads): + atmos_refract, numthreads, sst=False): """Calculate the solar position assuming unixtime is a numpy array. Note this function will not work if the solar position functions were compiled with numba. @@ -993,6 +1034,10 @@ def solar_position_numpy(unixtime, lat, lon, elev, pressure, temp, delta_t, v = apparent_sidereal_time(v0, delta_psi, epsilon) alpha = geocentric_sun_right_ascension(lamd, epsilon, beta) delta = geocentric_sun_declination(lamd, epsilon, beta) + if sst: + return v, alpha, delta + m = sun_mean_longitude(jme) + eot = equation_of_time(m, alpha, delta_psi, epsilon) H = local_hour_angle(v, lon, alpha) xi = equatorial_horizontal_parallax(R) u = uterm(lat) @@ -1010,19 +1055,19 @@ def solar_position_numpy(unixtime, lat, lon, elev, pressure, temp, delta_t, theta0 = topocentric_zenith_angle(e0) gamma = topocentric_astronomers_azimuth(H_prime, delta_prime, lat) phi = topocentric_azimuth_angle(gamma) - return theta, theta0, e, e0, phi + return theta, theta0, e, e0, phi, eot def solar_position(unixtime, lat, lon, elev, pressure, temp, delta_t, - atmos_refract, numthreads=8,): + atmos_refract, numthreads=8, sst=False): """ - Calculate the solar position using a translation of the - NREL SPA code. + Calculate the solar position using the + NREL SPA algorithm described in [1]. If numba is installed, the functions can be compiled and the code runs quickly. If not, the functions - still evaluate but use numpy instead + still evaluate but use numpy instead. Parameters ---------- @@ -1042,12 +1087,15 @@ def solar_position(unixtime, lat, lon, elev, pressure, temp, delta_t, temp : int or float avg. yearly temperature at location in degrees C; used for atmospheric correction - numthreads: int, optional - Number of threads to use for computation if numba>=0.17 - is installed. delta_t : float, optional Difference between terrestrial time and UT1. By default, use USNO historical data and predictions + atmos_refrac : float, optional + The approximate atmospheric refraction (in degrees) + at sunrise and sunset. + numthreads: int, optional + Number of threads to use for computation if numba>=0.17 + is installed. Returns ------- @@ -1056,7 +1104,8 @@ def solar_position(unixtime, lat, lon, elev, pressure, temp, delta_t, zenith, elevation, apparent_elevation, - azimuth + azimuth, + equation_of_time References ---------- @@ -1070,7 +1119,8 @@ def solar_position(unixtime, lat, lon, elev, pressure, temp, delta_t, result = do_calc(unixtime, lat, lon, elev, pressure, - temp, delta_t, atmos_refract, numthreads) + temp, delta_t, atmos_refract, numthreads, + sst) if not isinstance(result, np.ndarray): try: @@ -1079,3 +1129,111 @@ def solar_position(unixtime, lat, lon, elev, pressure, temp, delta_t, pass return result + + +def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads): + """ + Calculate the sun transit, sunrise, and sunset + for a set of dates at a given location. + + Parameters + ---------- + dates : array + Numpy array of ints/floats corresponding to the Unix time + for the dates of interest, must be midnight UTC (00:00+00:00) + on the day of interest. + lat : float + Latitude of location to perform calculation for + lon : float + Longitude of location + delta_t : float + Difference between terrestrial time and UT. USNO has tables. + numthreads : int + Number to threads to use for calculation (if using numba) + + Returns + ------- + tuple : (transit, sunrise, sunset) localized to UTC + + """ + + if ((dates % 86400) != 0.0).any(): + raise ValueError('Input dates must be at 00:00 UTC') + + utday = (dates // 86400) * 86400 + ttday0 = utday - delta_t + ttdayn1 = ttday0 - 86400 + ttdayp1 = ttday0 + 86400 + + #index 0 is v, 1 is alpha, 2 is delta + utday_res = solar_position(utday, 0, 0, 0, 0, 0, delta_t, + 0, numthreads, sst=True) + v = utday_res[0] + + ttday0_res = solar_position(ttday0, 0, 0, 0, 0, 0, delta_t, + 0, numthreads, sst=True) + ttdayn1_res = solar_position(ttdayn1, 0, 0, 0, 0, 0, delta_t, + 0, numthreads, sst=True) + ttdayp1_res = solar_position(ttdayp1, 0, 0, 0, 0, 0, delta_t, + 0, numthreads, sst=True) + m0 = (ttday0_res[1] - lon - v) / 360 + cos_arg = ((np.sin(np.radians(-0.8333)) - np.sin(np.radians(lat)) + * np.sin(np.radians(ttday0_res[2]))) / + (np.cos(np.radians(lat)) * np.cos(np.radians(ttday0_res[2])))) + cos_arg[abs(cos_arg) > 1] = np.nan + H0 = np.degrees(np.arccos(cos_arg)) % 180 + + m = np.empty((3, len(utday))) + m[0] = m0 % 1 + m[1] = (m[0] - H0 / 360) #% 1 + m[2] = (m[0] + H0 / 360) #% 1 + + # need to account for fractions of day that may be the next or previous + # day in UTC + add_a_day = m[2] >= 1 + sub_a_day = m[1] < 0 + m[1] = m[1] %1 + m[2] = m[2] %1 + vs = v + 360.985647 * m + n = m + delta_t / 86400 + + a = ttday0_res[1] - ttdayn1_res[1] + a[abs(a) > 2] = a[abs(a) > 2] % 1 + ap = ttday0_res[2] - ttdayn1_res[2] + ap[abs(ap) > 2] = ap[abs(ap) > 2] % 1 + b = ttdayp1_res[1] - ttday0_res[1] + b[abs(b) > 2] = b[abs(b) > 2] % 1 + bp = ttdayp1_res[2] - ttday0_res[2] + bp[abs(bp) > 2] = bp[abs(bp) > 2] % 1 + c = b - a + cp = bp - ap + + alpha_prime = ttday0_res[1] + (n * (a + b + c * n)) / 2 + delta_prime = ttday0_res[2] + (n * (ap + bp + cp * n)) / 2 + Hp = (vs + lon - alpha_prime) % 360 + Hp[Hp >= 180] = Hp[Hp >= 180] - 360 + + h = np.degrees(np.arcsin(np.sin(np.radians(lat)) * + np.sin(np.radians(delta_prime)) + + np.cos(np.radians(lat)) * + np.cos(np.radians(delta_prime)) + * np.cos(np.radians(Hp)))) + + T = (m[0] - Hp[0] / 360) * 86400 + R = (m[1] + (h[1] + 0.8333) / (360 * np.cos(np.radians(delta_prime[1])) * + np.cos(np.radians(lat)) * + np.sin(np.radians(Hp[1])))) * 86400 + S = (m[2] + (h[2] + 0.8333) / (360 * np.cos(np.radians(delta_prime[2])) * + np.cos(np.radians(lat)) * + np.sin(np.radians(Hp[2])))) * 86400 + + + S[add_a_day] += 86400 + R[sub_a_day] -= 86400 + + transit = T + utday + sunrise = R + utday + sunset = S + utday + + + return transit, sunrise, sunset diff --git a/pvlib/test/test_solarposition.py b/pvlib/test/test_solarposition.py index a2f087c5c7..8123f7bfed 100644 --- a/pvlib/test/test_solarposition.py +++ b/pvlib/test/test_solarposition.py @@ -4,6 +4,7 @@ import datetime import numpy as np +import numpy.testing as npt import pandas as pd from nose.tools import raises, assert_almost_equals @@ -32,7 +33,7 @@ def test_spa_physical(): times = pd.date_range(datetime.datetime(2003,10,17,12,30,30), periods=1, freq='D') try: - ephem_data = solarposition.spa(times, golden_mst, pressure=82000, + ephem_data = solarposition.spa_c(times, golden_mst, pressure=82000, temperature=11).ix[0] except ImportError: raise SkipTest @@ -43,9 +44,10 @@ def test_spa_physical(): def test_spa_physical_dst(): - times = pd.date_range(datetime.datetime(2003,10,17,13,30,30), periods=1, freq='D') + times = pd.date_range(datetime.datetime(2003,10,17,13,30,30), periods=1, + freq='D') try: - ephem_data = solarposition.spa(times, golden, pressure=82000, + ephem_data = solarposition.spa_c(times, golden, pressure=82000, temperature=11).ix[0] except ImportError: raise SkipTest @@ -57,7 +59,8 @@ def test_spa_physical_dst(): def test_spa_localization(): try: - assert_frame_equal(solarposition.spa(times, tus), solarposition.spa(times_localized, tus)) + assert_frame_equal(solarposition.spa_c(times, tus), + solarposition.spa_c(times_localized, tus)) except ImportError: raise SkipTest @@ -129,6 +132,42 @@ def test_spa_python_localization(): solarposition.spa_python(times_localized, tus)) +def test_get_sun_rise_set_transit(): + south = Location(-35.0, 0.0, tz='UTC') + times = pd.DatetimeIndex([datetime.datetime(1996, 7, 5, 0), + datetime.datetime(2004, 12, 4, 0)] + ).tz_localize('UTC') + sunrise = pd.DatetimeIndex([datetime.datetime(1996, 7, 5, 7, 8, 15, 471676), + datetime.datetime(2004, 12, 4, 4, 38, 57, 27416)] + ).tz_localize('UTC').tolist() + sunset = pd.DatetimeIndex([datetime.datetime(1996, 7, 5, 17, 1, 4, 479889), + datetime.datetime(2004, 12, 4, 19, 2, 2, 499704)] + ).tz_localize('UTC').tolist() + result = solarposition.get_sun_rise_set_transit(times, south, + delta_t=64.0) + frame = pd.DataFrame({'sunrise':sunrise, 'sunset':sunset}, index=times) + del result['transit'] + assert_frame_equal(frame, result) + + + # tests from USNO + # Golden + golden = Location(39.0, -105.0, tz='MST') + times = pd.DatetimeIndex([datetime.datetime(2015, 1, 2), + datetime.datetime(2015, 8, 2),] + ).tz_localize('MST') + sunrise = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 7, 19, 2, 225169), + datetime.datetime(2015, 8, 2, 5, 1, 26, 963145) + ]).tz_localize('MST').tolist() + sunset = pd.DatetimeIndex([datetime.datetime(2015, 1, 2, 16, 49, 10, 13145), + datetime.datetime(2015, 8, 2, 19, 11, 31, 816401) + ]).tz_localize('MST').tolist() + result = solarposition.get_sun_rise_set_transit(times, golden, delta_t=64.0) + frame = pd.DataFrame({'sunrise':sunrise, 'sunset':sunset}, index=times) + del result['transit'] + assert_frame_equal(frame, result) + + def test_pyephem_physical(): times = pd.date_range(datetime.datetime(2003,10,17,12,30,30), periods=1, freq='D') ephem_data = solarposition.pyephem(times, golden_mst, pressure=82000, temperature=11).ix[0] diff --git a/pvlib/test/test_spa.py b/pvlib/test/test_spa.py index 457c315268..0eb1652a67 100644 --- a/pvlib/test/test_spa.py +++ b/pvlib/test/test_spa.py @@ -1,4 +1,5 @@ import os +import datetime as dt import logging pvl_logger = logging.getLogger('pvlib') try: @@ -228,7 +229,88 @@ def test_solar_position(self): npt.assert_almost_equal( np.array([[theta, theta0, e, e0, Phi]]).T, self.spa.solar_position( unixtimes, lat, lon, elev, pressure, temp, delta_t, - atmos_refract), 5) + atmos_refract)[:-1], 5) + npt.assert_almost_equal( + np.array([[v, alpha, delta]]).T, self.spa.solar_position( + unixtimes, lat, lon, elev, pressure, temp, delta_t, + atmos_refract, sst=True)[:3], 5) + + def test_equation_of_time(self): + eot = 14.64 + M = self.spa.sun_mean_longitude(JME) + assert_almost_equals(eot, self.spa.equation_of_time( + M, alpha, dPsi, epsilon), 2) + + + def test_transit_sunrise_sunset(self): + # tests at greenwich + times = pd.DatetimeIndex([dt.datetime(1996, 7, 5, 0), + dt.datetime(2004, 12, 4, 0)] + ).tz_localize('UTC').astype(int)*1.0/10**9 + sunrise = pd.DatetimeIndex([dt.datetime(1996, 7, 5, 7, 8, 15), + dt.datetime(2004, 12, 4, 4, 38, 57)] + ).tz_localize('UTC').astype(int)*1.0/10**9 + sunset = pd.DatetimeIndex([dt.datetime(1996, 7, 5, 17, 1, 4), + dt.datetime(2004, 12, 4, 19, 2, 2)] + ).tz_localize('UTC').astype(int)*1.0/10**9 + result = self.spa.transit_sunrise_sunset(times, -35.0, 0.0, 64.0, 1) + npt.assert_almost_equal(sunrise/1e3, result[1]/1e3, 3) + npt.assert_almost_equal(sunset/1e3, result[2]/1e3, 3) + + + times = pd.DatetimeIndex([dt.datetime(1994, 1, 2),] + ).tz_localize('UTC').astype(int)*1.0/10**9 + sunset = pd.DatetimeIndex([dt.datetime(1994, 1, 2, 16, 59, 55),] + ).tz_localize('UTC').astype(int)*1.0/10**9 + sunrise = pd.DatetimeIndex([dt.datetime(1994, 1, 2, 7, 8, 12),] + ).tz_localize('UTC').astype(int)*1.0/10**9 + result = self.spa.transit_sunrise_sunset(times, 35.0, 0.0, 64.0, 1) + npt.assert_almost_equal(sunrise/1e3, result[1]/1e3, 3) + npt.assert_almost_equal(sunset/1e3, result[2]/1e3, 3) + + # tests from USNO + # Golden + times = pd.DatetimeIndex([dt.datetime(2015, 1, 2), + dt.datetime(2015, 4, 2), + dt.datetime(2015, 8, 2), + dt.datetime(2015, 12, 2),], + ).tz_localize('UTC').astype(int)*1.0/10**9 + sunrise = pd.DatetimeIndex([dt.datetime(2015, 1, 2, 7, 19), + dt.datetime(2015, 4, 2, 5, 43), + dt.datetime(2015, 8, 2, 5, 1), + dt.datetime(2015, 12, 2, 7, 1),], + ).tz_localize('MST').astype(int)*1.0/10**9 + sunset = pd.DatetimeIndex([dt.datetime(2015, 1, 2, 16, 49), + dt.datetime(2015, 4, 2, 18, 24), + dt.datetime(2015, 8, 2, 19, 10), + dt.datetime(2015, 12, 2, 16, 38),], + ).tz_localize('MST').astype(int)*1.0/10**9 + result = self.spa.transit_sunrise_sunset(times, 39.0, -105.0, 64.0, 1) + npt.assert_almost_equal(sunrise/1e3, result[1]/1e3, 1) + npt.assert_almost_equal(sunset/1e3, result[2]/1e3, 1) + + # Beijing + times = pd.DatetimeIndex([dt.datetime(2015, 1, 2), + dt.datetime(2015, 4, 2), + dt.datetime(2015, 8, 2), + dt.datetime(2015, 12, 2),], + ).tz_localize('UTC').astype(int)*1.0/10**9 + sunrise = pd.DatetimeIndex([dt.datetime(2015, 1, 2, 7, 36), + dt.datetime(2015, 4, 2, 5, 58), + dt.datetime(2015, 8, 2, 5, 13), + dt.datetime(2015, 12, 2, 7, 17),], + ).tz_localize('Asia/Shanghai' + ).astype(int)*1.0/10**9 + sunset = pd.DatetimeIndex([dt.datetime(2015, 1, 2, 17, 0), + dt.datetime(2015, 4, 2, 18, 39), + dt.datetime(2015, 8, 2, 19, 28), + dt.datetime(2015, 12, 2, 16, 50),], + ).tz_localize('Asia/Shanghai' + ).astype(int)*1.0/10**9 + result = self.spa.transit_sunrise_sunset(times, 39.917, 116.383, 64.0,1) + npt.assert_almost_equal(sunrise/1e3, result[1]/1e3, 1) + npt.assert_almost_equal(sunset/1e3, result[2]/1e3, 1) + class NumpySpaTest(unittest.TestCase, SpaBase): @@ -267,6 +349,16 @@ def tearDownClass(self): def test_julian_day(self): assert_almost_equals(JD, self.spa.julian_day(unixtimes[0]), 6) + def test_solar_position_singlethreaded(self): + npt.assert_almost_equal( + np.array([[theta, theta0, e, e0, Phi]]).T, self.spa.solar_position( + unixtimes, lat, lon, elev, pressure, temp, delta_t, + atmos_refract, numthreads=1)[:-1], 5) + npt.assert_almost_equal( + np.array([[v, alpha, delta]]).T, self.spa.solar_position( + unixtimes, lat, lon, elev, pressure, temp, delta_t, + atmos_refract, numthreads=1, sst=True)[:3], 5) + def test_solar_position_multithreaded(self): result = np.array([theta, theta0, e, e0, Phi]) nresult = np.array([result, result, result]).T @@ -275,5 +367,12 @@ def test_solar_position_multithreaded(self): nresult , self.spa.solar_position( times, lat, lon, elev, pressure, temp, delta_t, - atmos_refract, numthreads=8), 5) + atmos_refract, numthreads=8)[:-1], 5) + result = np.array([v, alpha, delta]) + nresult = np.array([result, result, result]).T + npt.assert_almost_equal( + nresult + , self.spa.solar_position( + times, lat, lon, elev, pressure, temp, delta_t, + atmos_refract, numthreads=8, sst=True)[:3], 5) diff --git a/setup.py b/setup.py index 8f157165bf..912882a085 100644 --- a/setup.py +++ b/setup.py @@ -35,7 +35,6 @@ 'pandas >= 0.13.1', 'pytz', 'six', - 'pyephem', ], 'scripts': [], 'include_package_data': True