add 'nrel' spa option for ET irradiance, fix doy - 1 bug (#215)

* use extraradiation function in disc. creates doy-1 offset problem. add invalid model to extrarad

* switch to doy-1. need new tests

* add numba-capable nrel_earthsun_distance

* add nrel option to extraradiation

* update irradiance notebook with new dni_extra method

* fix docstrings

* add check_less_precise=2 to failing comparisons

* remove unnecessary and slow date conversion

* fix tests

* fix bad rebase

* comment

* update whatsnew

Author: wholmgren

docs/sphinx/source/whatsnew/v0.4.0.txt

@@ -28,13 +28,21 @@ Enhancements
 * Adds the PVWatts DC, AC, and system losses model. (:issue:`195`)
 * Improve PEP8 conformity in irradiance module. (:issue:`214`)
 * irradiance.disc is up to 10x faster. (:issue:`214`)
+* Add solarposition.nrel_earthsun_distance function and option to
+  calculate extraterrestrial radiation using the NREL solar position
+  algorithm. (:issue:`211`, :issue:`215`)
 
 
 Bug fixes
 ~~~~~~~~~
 
 * dirint function yielded the wrong results for non-sea-level pressures.
   Fixed. (:issue:`212`)
+* Fixed a bug in the day angle calculation used by the 'spencer' and 'asce'
+  extraterrestrial radiation options. Most modeling results will be changed
+  by less than 1 part in 1000. (:issue:`211`)
+* irradiance.extraradiation now raises a ValueError for invalid method
+  input. It previously failed silently. (:issue:`215`)
 
 
 Documentation

docs/tutorials/irradiance.ipynb

@@ -35,32 +35,35 @@
                    'dirty steel': 0.08}
 
 
-def extraradiation(datetime_or_doy, solar_constant=1366.1, method='spencer'):
+def extraradiation(datetime_or_doy, solar_constant=1366.1, method='spencer',
+                   **kwargs):
     """
     Determine extraterrestrial radiation from day of year.
 
     Parameters
     ----------
     datetime_or_doy : int, float, array, pd.DatetimeIndex
-        Day of year, array of days of year e.g. pd.DatetimeIndex.dayofyear,
-        or pd.DatetimeIndex.
+        Day of year, array of days of year e.g.
+        pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
 
     solar_constant : float
         The solar constant.
 
     method : string
         The method by which the ET radiation should be calculated.
-        Options include ``'pyephem', 'spencer', 'asce'``.
+        Options include ``'pyephem', 'spencer', 'asce', 'nrel'``.
+
+    kwargs :
+        Passed to solarposition.nrel_earthsun_distance
 
     Returns
     -------
-    float or Series
-
+    dni_extra : float, array, or Series
         The extraterrestrial radiation present in watts per square meter
-        on a surface which is normal to the sun. Ea is of the same size as the
-        input doy.
+        on a surface which is normal to the sun. Ea is of the same size
+        as the input doy.
 
-        'pyephem' always returns a series.
+        'pyephem' and 'nrel' always return a Series.
 
     Notes
     -----
@@ -69,8 +72,8 @@ def extraradiation(datetime_or_doy, solar_constant=1366.1, method='spencer'):
 
     References
     ----------
-    [1] M. Reno, C. Hansen, and J. Stein, "Global Horizontal Irradiance Clear
-    Sky Models: Implementation and Analysis", Sandia National
+    [1] M. Reno, C. Hansen, and J. Stein, "Global Horizontal Irradiance
+    Clear Sky Models: Implementation and Analysis", Sandia National
     Laboratories, SAND2012-2389, 2012.
 
     [2] <http://solardat.uoregon.edu/SolarRadiationBasics.html>,
@@ -102,8 +105,9 @@ def extraradiation(datetime_or_doy, solar_constant=1366.1, method='spencer'):
         doy = datetime_or_doy
         input_to_datetimeindex = _array_to_datetimeindex
 
-    B = (2 * np.pi / 365) * doy
+    B = (2. * np.pi / 365.) * (doy - 1)
 
+    method = method.lower()
     if method == 'asce':
         pvl_logger.debug('Calculating ET rad using ASCE method')
         RoverR0sqrd = 1 + 0.033 * np.cos(B)
@@ -115,6 +119,12 @@ def extraradiation(datetime_or_doy, solar_constant=1366.1, method='spencer'):
         pvl_logger.debug('Calculating ET rad using pyephem method')
         times = input_to_datetimeindex(datetime_or_doy)
         RoverR0sqrd = solarposition.pyephem_earthsun_distance(times) ** (-2)
+    elif method == 'nrel':
+        times = input_to_datetimeindex(datetime_or_doy)
+        RoverR0sqrd = \
+            solarposition.nrel_earthsun_distance(times, **kwargs) ** (-2)
+    else:
+        raise ValueError('Invalid method: %s', method)
 
     Ea = solar_constant * RoverR0sqrd
 
@@ -1075,7 +1085,7 @@ def perez(surface_tilt, surface_azimuth, dhi, dni, dni_extra,
     return sky_diffuse
 
 
-def disc(ghi, zenith, times, pressure=101325):
+def disc(ghi, zenith, datetime_or_doy, pressure=101325):
     """
     Estimate Direct Normal Irradiance from Global Horizontal Irradiance
     using the DISC model.
@@ -1093,7 +1103,9 @@ def disc(ghi, zenith, times, pressure=101325):
         True (not refraction-corrected) solar zenith angles in decimal
         degrees.
 
-    times : DatetimeIndex
+    datetime_or_doy : int, float, array, pd.DatetimeIndex
+        Day of year or array of days of year e.g.
+        pd.DatetimeIndex.dayofyear, or pd.DatetimeIndex.
 
     pressure : numeric
         Site pressure in Pascal.
@@ -1128,18 +1140,8 @@ def disc(ghi, zenith, times, pressure=101325):
     dirint
     """
 
-    # in principle, the dni_extra calculation could be done by
-    # pvlib's function. However, this is the algorithm used in
-    # the DISC paper
-
-    doy = times.dayofyear
-
-    dayangle = 2. * np.pi*(doy - 1) / 365
-
-    re = (1.00011 + 0.034221*np.cos(dayangle) + 0.00128*np.sin(dayangle) +
-          0.000719*np.cos(2.*dayangle) + 7.7e-5*np.sin(2.*dayangle))
-
-    I0 = re * 1370.
+    # this is the I0 calculation from the reference
+    I0 = extraradiation(datetime_or_doy, 1370, 'spencer')
     I0h = I0 * np.cos(np.radians(zenith))
 
     am = atmosphere.relativeairmass(zenith, model='kasten1966')
@@ -1176,8 +1178,8 @@ def disc(ghi, zenith, times, pressure=101325):
     output['kt'] = kt
     output['airmass'] = am
 
-    if isinstance(times, pd.DatetimeIndex):
-        output = pd.DataFrame(output, index=times)
+    if isinstance(datetime_or_doy, pd.DatetimeIndex):
+        output = pd.DataFrame(output, index=datetime_or_doy)
 
     return output
 

pvlib/irradiance.py

@@ -771,3 +771,54 @@ 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 the
+    NREL SPA algorithm described in [1].
+
+    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.
+
+    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.
+    """
+    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, ])
+
+    unixtime = time.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/solarposition.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]
@@ -915,9 +916,12 @@ def solar_position_loop(unixtime, loc_args, out):
         jc = julian_century(jd)
         jce = julian_ephemeris_century(jde)
         jme = julian_ephemeris_millennium(jce)
+        R = heliocentric_radius_vector(jme)
+        if esd:
+            out[0, i] = R
+            continue
         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)
@@ -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.
@@ -1017,9 +1032,11 @@ def solar_position_numpy(unixtime, lat, lon, elev, pressure, temp, delta_t,
     jc = julian_century(jd)
     jce = julian_ephemeris_century(jde)
     jme = julian_ephemeris_millennium(jce)
+    R = heliocentric_radius_vector(jme)
+    if esd:
+        return (R, )
     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)
@@ -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,37 @@ def transit_sunrise_sunset(dates, lat, lon, delta_t, numthreads):
     sunset = S + utday
 
     return transit, sunrise, sunset
+
+
+def earthsun_distance(unixtime, delta_t, numthreads):
+    """
+    Calculates the distance from the earth to the sun using the
+    NREL SPA algorithm described in [1].
+
+    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.
+
+    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.
+    """
+
+    R = solar_position(unixtime, 0, 0, 0, 0, 0, delta_t,
+                       0, numthreads, esd=True)[0]
+
+    return R