Add a new and faster way to calculate the solarposition

.stickler.yml

@@ -2,7 +2,7 @@ linters:
     flake8:
         python: 3
         max-line-length: 79
-        ignore: E201,E241,E226
+        ignore: E201,E241,E226,W503,W504
 files:
     ignore:
         - 'pvlib/_version.py'

docs/sphinx/source/api.rst

@@ -45,6 +45,7 @@ unless you know that you need a different function.
    solarposition.ephemeris
    solarposition.pyephem
    solarposition.spa_c
+   solarposition.spencer
 
 
 Additional functions for quantities closely related to solar position.
@@ -55,6 +56,7 @@ Additional functions for quantities closely related to solar position.
    solarposition.calc_time
    solarposition.pyephem_earthsun_distance
    solarposition.nrel_earthsun_distance
+   solarposition.datetime2julian
    spa.calculate_deltat
 
 

docs/sphinx/source/whatsnew/v0.7.0.rst

@@ -8,8 +8,11 @@ recommend all users of v0.6.3 upgrade to this release.
 
 **Python 2.7 support ended on June 1, 2019**. (:issue:`501`)
 
+- Add a faster way to calculate the solar position (`spencer`)
+
 
 Contributors
 ~~~~~~~~~~~~
 * Mark Campanellli (:ghuser:`markcampanelli`)
 * Will Holmgren (:ghuser:`wholmgren`)
+* Daniel Lassahn (:ghuser:`meteoDaniel`)

pvlib/solarposition.py

@@ -7,6 +7,7 @@
 # Will Holmgren (@wholmgren), University of Arizona, 2014
 # Tony Lorenzo (@alorenzo175), University of Arizona, 2015
 # Cliff hansen (@cwhanse), Sandia National Laboratories, 2018
+# Daniel Lassahn (@meteoDaniel), meteocontrol Energy and Weather Services, 2019
 
 from __future__ import division
 import os
@@ -18,17 +19,16 @@
         from imp import reload
     except ImportError:
         pass
-
 import numpy as np
 import pandas as pd
 import warnings
 
 from pvlib import atmosphere
-from pvlib.tools import datetime_to_djd, djd_to_datetime
+from pvlib.tools import datetime_to_djd, djd_to_datetime, datetime_to_julian
 from pvlib._deprecation import deprecated
 
-
 NS_PER_HR = 1.e9 * 3600.  # nanoseconds per hour
+JULIAN_YEARS = 365.2425
 
 
 def get_solarposition(time, latitude, longitude,
@@ -67,6 +67,9 @@ def get_solarposition(time, latitude, longitude,
 
         'nrel_c' uses the NREL SPA C code [3]: :py:func:`spa_c`
 
+        'spencer' uses the  Spencer formula [4] :py:func:`spencer`
+
+
     temperature : float, default 12
         Degrees C.
 
@@ -81,6 +84,9 @@ def get_solarposition(time, latitude, longitude,
     solar radiation applications. Solar Energy, vol. 81, no. 6, p. 838, 2007.
 
     [3] NREL SPA code: http://rredc.nrel.gov/solar/codesandalgorithms/spa/
+
+    [4] Spencer (1972) can be found in "Solar energy fundamentals and
+    modeling techniques" from Zekai Sen
     """
 
     if altitude is None and pressure is None:
@@ -112,8 +118,10 @@ def get_solarposition(time, latitude, longitude,
                            pressure=pressure,
                            temperature=temperature, **kwargs)
     elif method == 'ephemeris':
-        ephem_df = ephemeris(time, latitude, longitude, pressure, temperature,
-                             **kwargs)
+        ephem_df = ephemeris(time, latitude, longitude, pressure, temperature)
+    elif method == 'spencer':
+        ephem_df = spencer(time, latitude, longitude)
+
     else:
         raise ValueError('Invalid solar position method')
 
@@ -495,7 +503,7 @@ def sun_rise_set_transit_ephem(times, latitude, longitude,
 
     Parameters
     ----------
-    time : pandas.DatetimeIndex
+    times : :class:`pandas.DatetimeIndex`
         Must be localized
     latitude : float
         Latitude in degrees, positive north of equator, negative to south
@@ -1445,3 +1453,80 @@ def sun_rise_set_transit_geometric(times, latitude, longitude, declination,
     sunset = _local_times_from_hours_since_midnight(times, sunset_hour)
     transit = _local_times_from_hours_since_midnight(times, transit_hour)
     return sunrise, sunset, transit
+
+
+def spencer(times, latitude, longitude):
+    """
+    Calculate the solar position using a formulation by Spencer 1971/1972.
+
+    Parameters
+    ----------
+    times : pandas.DatetimeIndex`
+        Corresponding timestamps, must be localized to the timezone for the
+        ``latitude`` and ``longitude``.
+    latitude : float
+        Latitude in degrees, positive north of equator, negative to south
+    longitude : float
+        Longitude in degrees, positive east of prime meridian, negative to west
+
+    Returns
+    -------
+    DataFrame
+        The DataFrame will have the following columns:
+        zenith (degrees),
+        elevation (degrees),
+        azimuth (degrees),
+        equation_of_time (seconds),
+        declination (degrees).
+
+    References
+    ----------
+    [1] Spencer (1972) can be found in
+    Sen, Zekai. Solar energy fundamentals and modeling techniques:
+     atmosphere, environment, climate change and renewable energy.
+     Springer Science & Business Media, 2008.
+    """
+
+    julians, julian_2000 = datetime_to_julian(times)
+    julians_2000 = np.asarray(julians, dtype=np.float) - julian_2000
+
+    latitude_radians = np.radians(latitude)
+    day_time = (julians_2000 % 1) * 24
+
+    declination = np.array(declination_spencer71(times.dayofyear))
+
+    # Equation of time (difference between standard time and solar time).
+    eot = np.array(equation_of_time_spencer71(times.dayofyear))
+
+    # True local time
+    tlt = (day_time + longitude / 15 + eot / 60) % 24 - 12
+
+    # Solar hour angle
+    ha = np.radians(tlt * 15)
+
+    # Calculate sun elevation.
+    sin_sun_elevation = (
+            np.sin(declination) * np.sin(latitude_radians) +
+            np.cos(declination) * np.cos(latitude_radians) * np.cos(ha)
+    )
+
+    # Compute the sun's elevation and zenith angle.
+    elevation = np.arcsin(sin_sun_elevation)
+    zenith = np.pi / 2 - elevation
+
+    # Compute the sun's azimuth angle.
+    y = -(np.sin(latitude_radians) * np.sin(elevation) - np.sin(declination)) \
+        / (np.cos(latitude_radians) * np.cos(elevation))
+    azimuth = np.arccos(y)
+
+    # Convert azimuth angle from 0-pi to 0-2pi.
+    tlt_filter = 0 <= tlt
+    azimuth[tlt_filter] = 2 * np.pi - azimuth[tlt_filter]
+
+    result = pd.DataFrame({'zenith': np.degrees(zenith),
+                           'elevation': np.degrees(elevation),
+                           'azimuth': np.degrees(azimuth),
+                           'declination': declination,
+                           'equation_of_time': eot},
+                          index=times)
+    return result

pvlib/spa.py

@@ -1161,7 +1161,7 @@ def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads):
 
     Parameters
     ----------
-    dates : array
+    dates : np.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.

pvlib/test/test_solarposition.py

@@ -754,8 +754,28 @@ def test_sun_rise_set_transit_geometric(expected_rise_set_spa, golden_mst):
                        atol=np.abs(expected_transit_error).max())
 
 
-# put numba tests at end of file to minimize reloading
+def test_get_solar_position_spencer(golden_mst):
+    """ test for the calculation based on spencer 1972 """
+    times = pd.date_range(datetime.datetime(2003, 10, 17, 12, 30, 30),
+                          periods=1, freq='D')
+    spencer_data = solarposition.spencer(times,
+                                         golden_mst.latitude,
+                                         golden_mst.longitude)
+
+    expected_solpos = golden_mst.get_solarposition(times)
 
+    np.testing.assert_array_almost_equal(spencer_data.zenith.values,
+                                         expected_solpos.zenith.values,
+                                         decimal=0)
+    np.testing.assert_array_almost_equal(spencer_data.azimuth.values,
+                                         expected_solpos.azimuth.values,
+                                         decimal=0)
+    np.testing.assert_array_almost_equal(spencer_data.elevation.values,
+                                         expected_solpos.elevation.values,
+                                         decimal=0)
+
+
+# put numba tests at end of file to minimize reloading
 @requires_numba
 def test_spa_python_numba_physical(expected_solpos, golden_mst):
     times = pd.date_range(datetime.datetime(2003, 10, 17, 12, 30, 30),

pvlib/test/test_tools.py

@@ -1,7 +1,14 @@
+import datetime
+
+import numpy as np
+import pandas as pd
 import pytest
 
 from pvlib import tools
 
+times = pd.date_range(start=datetime.datetime(2014, 6, 24),
+                      end=datetime.datetime(2014, 6, 26), freq='15Min')
+
 
 @pytest.mark.parametrize('keys, input_dict, expected', [
     (['a', 'b'], {'a': 1, 'b': 2, 'c': 3}, {'a': 1, 'b': 2}),
@@ -12,3 +19,21 @@
 def test_build_kwargs(keys, input_dict, expected):
     kwargs = tools._build_kwargs(keys, input_dict)
     assert kwargs == expected
+
+
+def test_datetime_to_julian():
+    """ test transformation from datetime to julians """
+    julians = tools.datetime_to_julian(pd.to_datetime(times))
+    np.testing.assert_array_almost_equal(np.array(julians[:10]),
+                                         np.array([
+                                             2456832.5,
+                                             2456832.5104166665,
+                                             2456832.5208333335,
+                                             2456832.53125,
+                                             2456832.5416666665,
+                                             2456832.5520833335,
+                                             2456832.5625,
+                                             2456832.5729166665,
+                                             2456832.5833333335,
+                                             2456832.59375])
+                                         )

pvlib/tools.py

@@ -9,6 +9,10 @@
 import pandas as pd
 import pytz
 
+JULIAN_2000 = 2451545
+DAY_SECONDS = 60 * 60 * 24
+DT_2000 = pd.to_datetime(dt.datetime(2000, 1, 1)).tz_localize('UTC')
+
 
 def cosd(angle):
     """
@@ -425,3 +429,30 @@ def _golden_sect_DataFrame(params, VL, VH, func):
             raise Exception("EXCEPTION:iterations exceeded maximum (50)")
 
     return func(df, 'V1'), df['V1']
+
+
+def datetime_to_julian(times):
+    """
+        Transforms pd.DateTimeIndex to Julian days
+
+    Parameters
+    ----------
+    times : pandas.DatetimeIndex
+        Corresponding timestamps, must be localized to the timezone for the
+        ``latitude`` and ``longitude``
+    Returns
+    -------
+    Float64Index
+        The float index contains julian dates
+    Float
+        julian 2000 in UTC referenced to local time
+    """
+    dt_2000 = DT_2000.tz_convert(times.tzinfo)
+    hours_difference = (DT_2000.tz_localize(None) -
+                        dt_2000.tz_localize(None)).seconds / 3600
+    julian_2000 = JULIAN_2000 - (hours_difference/24)
+    delta = times - dt_2000
+    delta_julians = (delta.seconds + delta.microseconds / 1e6)
+    return (
+        julian_2000 + delta.days + delta_julians / DAY_SECONDS
+    ), julian_2000