Add a new and faster way to calculate the solarposition

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_mc
 
 
 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.1.rst

@@ -0,0 +1,11 @@
+.. _whatsnew_0701:
+
+v0.7.1 (MONTH DAY, YEAR)
+---------------------
+
+- Add a faster way to calculate the solar position (`spencer_mc`)
+
+
+Contributors
+~~~~~~~~~~~~
+* 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,7 +19,6 @@
         from imp import reload
     except ImportError:
         pass
-
 import numpy as np
 import pandas as pd
 import warnings
@@ -27,8 +27,12 @@
 from pvlib.tools import datetime_to_djd, djd_to_datetime
 from pvlib._deprecation import deprecated
 
-
 NS_PER_HR = 1.e9 * 3600.  # nanoseconds per hour
+JULIAN_2000 = 2451544.5
+DT_2000 = dt.datetime(2000, 1, 1)
+TIMESTAMP_2000 = DT_2000.timestamp()
+DAY_SECONDS = 60 * 60 * 24
+JULIAN_YEARS = 365.2425
 
 
 def get_solarposition(time, latitude, longitude,
@@ -67,6 +71,9 @@ def get_solarposition(time, latitude, longitude,
 
         'nrel_c' uses the NREL SPA C code [3]: :py:func:`spa_c`
 
+        'spencer_mc' uses the  Spencer formula [4] :py:func:`spencer_mc`
+
+
     temperature : float, default 12
         Degrees C.
 
@@ -81,6 +88,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 +122,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_mc':
+        ephem_df = spencer_mc(time, latitude, longitude)
+
     else:
         raise ValueError('Invalid solar position method')
 
@@ -495,7 +507,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 +1457,118 @@ 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_mc(times, latitude, longitude):
+    """
+    Calculate the solar position using a python implementation of the
+    Spencer (1972) formulation provided by meteocontrol
+
+    Parameters
+    ----------
+    times : :class:`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),
+        eccentricity,
+        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 = datetime2julian(times)
+    julians_2000 = np.asarray(julians, dtype=np.float) - JULIAN_2000
+
+    lat = np.radians(latitude)
+
+    # Compute fractional year (gamma) in radians
+    gamma = 2 * np.pi * (julians_2000 % JULIAN_YEARS) / JULIAN_YEARS
+    cos_gamma = np.cos(gamma), np.cos(gamma * 2), np.cos(gamma * 3)
+    sin_gamma = np.sin(gamma), np.sin(gamma * 2), np.sin(gamma * 3)
+    day_time = (julians_2000 % 1) * 24
+
+    # Eccentricity: correction factor of the earth's orbit.
+    eccentricity = (
+            1.00011 +
+            0.034221 * cos_gamma[0] +
+            0.001280 * sin_gamma[0] +
+            0.000719 * cos_gamma[1] +
+            0.000077 * sin_gamma[1]
+    )
+
+    declination = np.array(declination_spencer71(times.dayofyear))
+
+    # Equation of time (difference between standard time and solar time).
+    eot = (0.000075 + 0.001868 * cos_gamma[0] - 0.032077 * sin_gamma[0]
+           - 0.014615 * cos_gamma[1] - 0.040849 * sin_gamma[1]) * 229.18
+
+    # 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(lat) +
+                         np.cos(declination) * np.cos(lat) * 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(lat) * np.sin(elevation) - np.sin(declination)) \
+        / (np.cos(lat) * 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),
+                           'eccentricity': eccentricity,
+                           'declination': declination,
+                           'equation_of_time': eot},
+                          index=times)
+    return result
+
+
+def datetime2julian(times):
+    """
+        Transforms pd.DateTimeIndex to Julian days
+
+    Parameters
+    ----------
+    times : :class:`pandas.DatetimeIndex`
+        Corresponding timestamps, must be localized to the timezone for the
+        ``latitude`` and ``longitude``
+    Returns
+    -------
+    Float64Index
+        The float index contains julian dates
+    """
+
+    delta = times - DT_2000
+    return (
+        JULIAN_2000 + delta.days +
+        (delta.seconds + delta.microseconds / 1e6) / DAY_SECONDS
+    )

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

@@ -753,7 +753,78 @@ def test_sun_rise_set_transit_geometric(expected_rise_set_spa, golden_mst):
     assert np.allclose(test_transit, expected_transit,
                        atol=np.abs(expected_transit_error).max())
 
-
+
+def test_spencer_mc():
+    """ test for the calculation based on spencer 1972 """
+    latitude = 48.367073
+    longitude = 10.868378
+    assert_frame_equal(solarposition.spencer_mc(times[:10],
+                                                latitude,
+                                                longitude),
+                       pd.DataFrame([[107.60190741, -17.60190741, 9.95454033,
+                                      0.96706069, 0.40900622, -2.0908958],
+                                     [107.09533954, -17.09533954, 13.51934494,
+                                      0.96705971,
+                                      0.4090031, -2.09315672],
+                                     [106.43889684, -16.43889684, 17.04622292,
+                                      0.96705874,
+                                      0.40899997, -2.0954175],
+                                     [105.63682427, -15.63682427, 20.52704485,
+                                      0.96705777,
+                                      0.40899683, -2.09767815],
+                                     [104.69411613, -14.69411613, 23.95490053,
+                                      0.9670568,
+                                      0.40899367, -2.09993867],
+                                     [103.61638656, -13.61638656, 27.32420466,
+                                      0.96705583,
+                                      0.4089905, -2.10219905],
+                                     [102.40973635, -12.40973635, 30.63074363,
+                                      0.96705486,
+                                      0.40898732, -2.1044593],
+                                     [101.08062211, -11.08062211, 33.87167089,
+                                      0.9670539,
+                                      0.40898412, -2.10671942],
+                                     [99.63573414, -9.63573414, 37.04545803,
+                                      0.96705293,
+                                      0.40898091, -2.10897939],
+                                     [98.0818867, -8.0818867, 40.15181379,
+                                      0.96705197,
+                                      0.40897769, -2.11123924]],
+                                    index=times,
+                                    columns=['zenith',
+                                             'elevation',
+                                             'azimuth',
+                                             'eccentricity',
+                                             'declination',
+                                             'equation_of_time'])
+                       )
+
+def test_datetime2julians():
+    """ test transformation from datetime to julians """
+    julians = solarposition.datetime2julian(pd.to_datetime(times))
+    np.testing.assert_array_almost_equal(julians,
+                                         np.array([2458119.5, 2458119.54166667,
+                                                   2458119.58333333,
+                                                   2458119.625,
+                                                   2458119.66666667,
+                                                   2458119.70833333,
+                                                   2458119.75, 2458119.79166667,
+                                                   2458119.83333333,
+                                                   2458119.875,
+                                                   2458119.91666667,
+                                                   2458119.95833333,
+                                                   2458120., 2458120.04166667,
+                                                   2458120.08333333,
+                                                   2458120.125,
+                                                   2458120.16666667,
+                                                   2458120.20833333,
+                                                   2458120.25, 2458120.29166667,
+                                                   2458120.33333333,
+                                                   2458120.375,
+                                                   2458120.41666667,
+                                                   2458120.45833333])
+                                         )
+
 # put numba tests at end of file to minimize reloading
 
 @requires_numba