add numba-capable nrel_earthsun_distance

Author: wholmgren

pvlib/solarposition.py

@@ -771,3 +771,49 @@ def pyephem_earthsun_distance(time):
         earthsun.append(sun.earth_distance)
 
     return pd.Series(earthsun, index=time)
+
+
+def nrel_earthsun_distance(time, how='numpy', delta_t=None, numthreads=4):
+    """
+    Calculates the distance from the earth to the sun using pyephem.
+
+    Parameters
+    ----------
+    time : pd.DatetimeIndex
+
+    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.
+
+    delta_t : float, optional
+        Difference between terrestrial time and UT1.
+        By default, use USNO historical data and predictions
+
+    numthreads : int, optional
+        Number of threads to use if how == 'numba'.
+
+    Returns
+    -------
+    R : pd.Series
+        Earth-sun distance in AU.
+    """
+    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(np.int64)/10**9
+
+    spa = _spa_python_import(how)
+
+    R = spa.earthsun_distance(unixtime, delta_t, numthreads)
+
+    R = pd.Series(R, index=time)
+
+    return R

pvlib/spa.py

@@ -907,6 +907,7 @@ def solar_position_loop(unixtime, loc_args, out):
     delta_t = loc_args[5]
     atmos_refract = loc_args[6]
     sst = loc_args[7]
+    esd = loc_args[8]
 
     for i in range(unixtime.shape[0]):
         utime = unixtime[i]
@@ -918,6 +919,9 @@ def solar_position_loop(unixtime, loc_args, out):
         L = heliocentric_longitude(jme)
         B = heliocentric_latitude(jme)
         R = heliocentric_radius_vector(jme)
+        if esd:
+            out[0, i] = R
+            continue
         Theta = geocentric_longitude(L)
         beta = geocentric_latitude(B)
         x0 = mean_elongation(jce)
@@ -970,14 +974,24 @@ def solar_position_loop(unixtime, loc_args, out):
 
 
 def solar_position_numba(unixtime, lat, lon, elev, pressure, temp, delta_t,
-                         atmos_refract, numthreads, sst=False):
+                         atmos_refract, numthreads, sst=False, esd=False):
     """Calculate the solar position using the numba compiled functions
     and multiple threads. Very slow if functions are not numba compiled.
     """
+    # these args are the same for each thread
     loc_args = np.array([lat, lon, elev, pressure, temp, delta_t,
-                         atmos_refract, sst])
+                         atmos_refract, sst, esd])
+
+    # construct dims x ulength array to put the results in
     ulength = unixtime.shape[0]
-    result = np.empty((6, ulength), dtype=np.float64)
+    if sst:
+        dims = 3
+    elif esd:
+        dims = 1
+    else:
+        dims = 6
+    result = np.empty((dims, ulength), dtype=np.float64)
+
     if unixtime.dtype != np.float64:
         unixtime = unixtime.astype(np.float64)
 
@@ -992,6 +1006,7 @@ def solar_position_numba(unixtime, lat, lon, elev, pressure, temp, delta_t,
         solar_position_loop(unixtime, loc_args, result)
         return result
 
+    # split the input and output arrays into numthreads chunks
     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)]
@@ -1006,7 +1021,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, sst=False):
+                         atmos_refract, numthreads, sst=False, esd=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.
@@ -1020,6 +1035,8 @@ def solar_position_numpy(unixtime, lat, lon, elev, pressure, temp, delta_t,
     L = heliocentric_longitude(jme)
     B = heliocentric_latitude(jme)
     R = heliocentric_radius_vector(jme)
+    if esd:
+        return (R, )
     Theta = geocentric_longitude(L)
     beta = geocentric_latitude(B)
     x0 = mean_elongation(jce)
@@ -1063,7 +1080,7 @@ def solar_position_numpy(unixtime, lat, lon, elev, pressure, temp, delta_t,
 
 
 def solar_position(unixtime, lat, lon, elev, pressure, temp, delta_t,
-                   atmos_refract, numthreads=8, sst=False):
+                   atmos_refract, numthreads=8, sst=False, esd=False):
 
     """
     Calculate the solar position using the
@@ -1100,6 +1117,11 @@ def solar_position(unixtime, lat, lon, elev, pressure, temp, delta_t,
     numthreads: int, optional
         Number of threads to use for computation if numba>=0.17
         is installed.
+    sst : bool
+        If True, return only data needed for sunrise, sunset, and transit
+        calculations.
+    esd : bool
+        If True, return only Earth-Sun distance in AU
 
     Returns
     -------
@@ -1126,7 +1148,7 @@ 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,
-                     sst)
+                     sst, esd)
 
     if not isinstance(result, np.ndarray):
         try:
@@ -1241,3 +1263,31 @@ def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads):
     sunset = S + utday
 
     return transit, sunrise, sunset
+
+
+def earthsun_distance(unixtime, delta_t, numthreads):
+    """
+    Calculate the sun transit, sunrise, and sunset
+    for a set of dates at a given location.
+
+    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(np.int64)/10**9
+    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
+    -------
+    R : array
+        Earth-Sun distance in AU.
+    """
+
+    R = solar_position(unixtime, 0, 0, 0, 0, 0, delta_t,
+                       0, numthreads, esd=True)[0]
+
+    return R

pvlib/test/test_solarposition.py

@@ -3,7 +3,8 @@
 import numpy as np
 import pandas as pd
 
-from pandas.util.testing import assert_frame_equal, assert_index_equal
+from pandas.util.testing import (assert_frame_equal, assert_series_equal,
+                                 assert_index_equal)
 from numpy.testing import assert_allclose
 import pytest
 
@@ -320,3 +321,19 @@ def test_get_solarposition_no_kwargs(expected_solpos):
     expected_solpos = np.round(expected_solpos, 2)
     ephem_data = np.round(ephem_data, 2)
     assert_frame_equal(expected_solpos, ephem_data[expected_solpos.columns])
+
+
+def test_nrel_earthsun_distance():
+    times = pd.DatetimeIndex([datetime.datetime(2015, 1, 2),
+                              datetime.datetime(2015, 8, 2),]
+                             ).tz_localize('MST')
+    result = solarposition.nrel_earthsun_distance(times, delta_t=64.0)
+    expected = pd.Series(np.array([0.983293144266, 1.01489789116]),
+                         index=times)
+    assert_series_equal(expected, result)
+
+    times = datetime.datetime(2015, 1, 2)
+    result = solarposition.nrel_earthsun_distance(times, delta_t=64.0)
+    expected = pd.Series(np.array([0.983293144266]),
+                         index=pd.DatetimeIndex([times, ]))
+    assert_series_equal(expected, result)

pvlib/test/test_spa.py

@@ -240,7 +240,6 @@ def test_equation_of_time(self):
         assert_almost_equal(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),
@@ -310,7 +309,12 @@ def test_transit_sunrise_sunset(self):
         assert_almost_equal(sunrise/1e3, result[1]/1e3, 1)
         assert_almost_equal(sunset/1e3, result[2]/1e3, 1)
 
-
+    def test_earthsun_distance(self):
+        times = (pd.date_range('2003-10-17 12:30:30', periods=1, freq='D')
+           .tz_localize('MST'))
+        unixtimes = times.tz_convert('UTC').astype(np.int64)*1.0/10**9
+        result = self.spa.earthsun_distance(unixtimes, 64.0, 1)
+        assert_almost_equal(R, result, 6)
 
 class NumpySpaTest(unittest.TestCase, SpaBase):
     """Import spa without compiling to numba then run tests"""