Issue 165: Added calculate_deltaT function into spa.py

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

@@ -16,15 +16,17 @@ API Changes
 ~~~~~~~~~~~
 
 * The run_model method of the ModelChain will use the weather parameter of all weather data instead of splitting it to irradiation and weather. The irradiation parameter still works but will be removed soon.(:issue:`239`)
+* delta_t kwarg is now 67.0 instead of None. IMPORTANT: Setting delta_t as None will break the code for Numba calculation. This will be fixed in a future version. (:issue:`165`)
 
 
 Enhancements
 ~~~~~~~~~~~~
 
 * Adding a complete_irradiance method to the ModelChain to make it possible to calculate missing irradiation data from the existing columns [beta] (:issue:`239`)
-
+* Added calculate_deltat method to the spa module to calculate the time difference between terrestrial time and UT1. Specifying a scalar is sufficient for most calculations. (:issue:`165`)
 
 Code Contributors
 ~~~~~~~~~~~~~~~~~
 
 * Uwe Krien
+* Volker Beutner

pvlib/solarposition.py

@@ -238,7 +238,7 @@ def _spa_python_import(how):
 
 
 def spa_python(time, latitude, longitude,
-               altitude=0, pressure=101325, temperature=12, delta_t=None,
+               altitude=0, pressure=101325, temperature=12, delta_t=67.0,
                atmos_refract=None, how='numpy', numthreads=4, **kwargs):
     """
     Calculate the solar position using a python implementation of the
@@ -261,7 +261,12 @@ def spa_python(time, latitude, longitude,
     temperature : int or float, optional
         avg. yearly air temperature in degrees C.
     delta_t : float, optional
+        If delta_t is None, uses spa.calculate_deltat
+        using time.year and time.month from pandas.DatetimeIndex.
+        For most simulations specifing delta_t is sufficient.
         Difference between terrestrial time and UT1.
+        *Note: delta_t = None will break code using nrel_numba,
+        this will be fixed in a future version.
         The USNO has historical and forecasted delta_t [3].
     atmos_refrac : float, optional
         The approximate atmospheric refraction (in degrees)
@@ -308,7 +313,7 @@ def spa_python(time, latitude, longitude,
     lon = longitude
     elev = 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):
@@ -321,6 +326,8 @@ def spa_python(time, latitude, longitude,
 
     spa = _spa_python_import(how)
 
+    delta_t = delta_t or spa.calculate_deltat(time.year, time.month)
+
     app_zenith, zenith, app_elevation, elevation, azimuth, eot = spa.solar_position(
         unixtime, lat, lon, elev, pressure, temperature, delta_t,
         atmos_refract, numthreads)
@@ -335,7 +342,7 @@ def spa_python(time, latitude, longitude,
 
 
 def get_sun_rise_set_transit(time, latitude, longitude, how='numpy',
-                             delta_t=None,
+                             delta_t=67.0,
                              numthreads=4):
     """
     Calculate the sunrise, sunset, and sun transit times using the
@@ -353,7 +360,12 @@ def get_sun_rise_set_transit(time, latitude, longitude, how='numpy',
     latitude : float
     longitude : float
     delta_t : float, optional
+        If delta_t is None, uses spa.calculate_deltat
+        using time.year and time.month from pandas.DatetimeIndex.
+        For most simulations specifing delta_t is sufficient.
         Difference between terrestrial time and UT1.
+        *Note: delta_t = None will break code using nrel_numba,
+        this will be fixed in a future version.
         By default, use USNO historical data and predictions
     how : str, optional
         Options are 'numpy' or 'numba'. If numba >= 0.17.0
@@ -380,7 +392,6 @@ def get_sun_rise_set_transit(time, latitude, longitude, how='numpy',
 
     lat = latitude
     lon = longitude
-    delta_t = delta_t or 67.0
 
     if not isinstance(time, pd.DatetimeIndex):
         try:
@@ -394,6 +405,8 @@ def get_sun_rise_set_transit(time, latitude, longitude, how='numpy',
 
     spa = _spa_python_import(how)
 
+    delta_t = delta_t or spa.calculate_deltat(time.year, time.month)
+
     transit, sunrise, sunset = spa.transit_sunrise_sunset(
         unixtime, lat, lon, delta_t, numthreads)
 
@@ -773,7 +786,7 @@ def pyephem_earthsun_distance(time):
     return pd.Series(earthsun, index=time)
 
 
-def nrel_earthsun_distance(time, how='numpy', delta_t=None, numthreads=4):
+def nrel_earthsun_distance(time, how='numpy', delta_t=67.0, numthreads=4):
     """
     Calculates the distance from the earth to the sun using the
     NREL SPA algorithm described in [1].
@@ -788,7 +801,12 @@ def nrel_earthsun_distance(time, how='numpy', delta_t=None, numthreads=4):
         to machine code and run them multithreaded.
 
     delta_t : float, optional
+        If delta_t is None, uses spa.calculate_deltat
+        using time.year and time.month from pandas.DatetimeIndex.
+        For most simulations specifing delta_t is sufficient.
         Difference between terrestrial time and UT1.
+        *Note: delta_t = None will break code using nrel_numba,
+        this will be fixed in a future version.
         By default, use USNO historical data and predictions
 
     numthreads : int, optional
@@ -805,7 +823,6 @@ def nrel_earthsun_distance(time, how='numpy', delta_t=None, numthreads=4):
     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:
@@ -817,6 +834,8 @@ def nrel_earthsun_distance(time, how='numpy', delta_t=None, numthreads=4):
 
     spa = _spa_python_import(how)
 
+    delta_t = delta_t or spa.calculate_deltat(time.year, time.month)
+
     R = spa.earthsun_distance(unixtime, delta_t, numthreads)
 
     R = pd.Series(R, index=time)

pvlib/spa.py

@@ -1109,7 +1109,12 @@ def solar_position(unixtime, lat, lon, elev, pressure, temp, delta_t,
         avg. yearly temperature at location in
         degrees C; used for atmospheric correction
     delta_t : float, optional
+        If delta_t is None, uses spa.calculate_deltat
+        using time.year and time.month from pandas.DatetimeIndex.
+        For most simulations specifing delta_t is sufficient.
         Difference between terrestrial time and UT1.
+        *Note: delta_t = None will break code using nrel_numba,
+        this will be fixed in a future version.
         By default, use USNO historical data and predictions
     atmos_refrac : float, optional
         The approximate atmospheric refraction (in degrees)
@@ -1297,3 +1302,134 @@ def earthsun_distance(unixtime, delta_t, numthreads):
                        0, numthreads, esd=True)[0]
 
     return R
+
+
+def calculate_deltat(year, month):
+    """Calculate the difference between Terrestrial Dynamical Time (TD)
+    and Universal Time (UT).
+
+    Note: This function is not yet compatible for calculations using
+    Numba.
+
+    Equations taken from http://eclipse.gsfc.nasa.gov/SEcat5/deltatpoly.html
+    """
+
+    plw = ' Deltat is unknown for years before -1999 and after 3000.'\
+          + ' Delta values will be calculated, but the calculations'\
+          + ' are not intended to be used for these years.'
+
+    try:
+        if np.any((year > 3000) | (year < -1999)):
+            pvl_logger.warning(plw)
+    except ValueError:
+        if (year > 3000) | (year < -1999):
+            pvl_logger.warning(plw)
+    except TypeError:
+        return 0
+
+    y = year + (month - 0.5)/12
+
+    deltat = np.where(year < -500,
+
+                      -20+32*((y-1820)/100)**2, 0)
+
+    deltat = np.where((-500 <= year) & (year < 500),
+
+                      10583.6-1014.41*(y/100)
+                      + 33.78311*(y/100)**2
+                      - 5.952053*(y/100)**3
+                      - 0.1798452*(y/100)**4
+                      + 0.022174192*(y/100)**5
+                      + 0.0090316521*(y/100)**6, deltat)
+
+    deltat = np.where((500 <= year) & (year < 1600),
+
+                      1574.2-556.01*((y-1000)/100)
+                      + 71.23472*((y-1000)/100)**2
+                      + 0.319781*((y-1000)/100)**3
+                      - 0.8503463*((y-1000)/100)**4
+                      - 0.005050998*((y-1000)/100)**5
+                      + 0.0083572073*((y-1000)/100)**6, deltat)
+
+    deltat = np.where((1600 <= year) & (year < 1700),
+
+                      120-0.9808*(y-1600)
+                      - 0.01532*(y-1600)**2
+                      + (y-1600)**3/7129, deltat)
+
+    deltat = np.where((1700 <= year) & (year < 1800),
+
+                      8.83+0.1603*(y-1700)
+                      - 0.0059285*(y-1700)**2
+                      + 0.00013336*(y-1700)**3
+                      - (y-1700)**4/1174000, deltat)
+
+    deltat = np.where((1800 <= year) & (year < 1860),
+
+                      13.72-0.332447*(y-1800)
+                      + 0.0068612*(y-1800)**2
+                      + 0.0041116*(y-1800)**3
+                      - 0.00037436*(y-1800)**4
+                      + 0.0000121272*(y-1800)**5
+                      - 0.0000001699*(y-1800)**6
+                      + 0.000000000875*(y-1800)**7, deltat)
+
+    deltat = np.where((1860 <= year) & (year < 1900),
+
+                      7.6+0.5737*(y-1860)
+                      - 0.251754*(y-1860)**2
+                      + 0.01680668*(y-1860)**3
+                      - 0.0004473624*(y-1860)**4
+                      + (y-1860)**5/233174, deltat)
+
+    deltat = np.where((1900 <= year) & (year < 1920),
+
+                      -2.79+1.494119*(y-1900)
+                      - 0.0598939*(y-1900)**2
+                      + 0.0061966*(y-1900)**3
+                      - 0.000197*(y-1900)**4, deltat)
+
+    deltat = np.where((1920 <= year) & (year < 1941),
+
+                      21.20+0.84493*(y-1920)
+                      - 0.076100*(y-1920)**2
+                      + 0.0020936*(y-1920)**3, deltat)
+
+    deltat = np.where((1941 <= year) & (year < 1961),
+
+                      29.07+0.407*(y-1950)
+                      - (y-1950)**2/233
+                      + (y-1950)**3/2547, deltat)
+
+    deltat = np.where((1961 <= year) & (year < 1986),
+
+                      45.45+1.067*(y-1975)
+                      - (y-1975)**2/260
+                      - (y-1975)**3/718, deltat)
+
+    deltat = np.where((1986 <= year) & (year < 2005),
+
+                      63.86+0.3345*(y-2000)
+                      - 0.060374*(y-2000)**2
+                      + 0.0017275*(y-2000)**3
+                      + 0.000651814*(y-2000)**4
+                      + 0.00002373599*(y-2000)**5, deltat)
+
+    deltat = np.where((2005 <= year) & (year < 2050),
+
+                      62.92+0.32217*(y-2000)
+                      + 0.005589*(y-2000)**2, deltat)
+
+    deltat = np.where((2050 <= year) & (year < 2150),
+
+                      -20+32*((y-1820)/100)**2
+                      - 0.5628*(2150-y), deltat)
+
+    deltat = np.where(year > 2150,
+
+                      -20+32*((y-1820)/100)**2, deltat)
+
+    deltat = np.asscalar(deltat) if np.isscalar(year) & np.isscalar(month)\
+        else deltat
+
+    return deltat

pvlib/test/test_irradiance.py

@@ -98,7 +98,7 @@ def test_grounddiffuse_albedo_invalid_surface():
 
 def test_grounddiffuse_albedo_surface():
     result = irradiance.grounddiffuse(40, ghi, surface_type='sand')
-    assert_allclose(result, [0, 3.731058, 48.778813, 12.035025])
+    assert_allclose(result, [0, 3.731058, 48.778813, 12.035025], atol=1e-4)
 
 
 def test_isotropic_float():
@@ -108,7 +108,7 @@ def test_isotropic_float():
 
 def test_isotropic_series():
     result = irradiance.isotropic(40, irrad_data['dhi'])
-    assert_allclose(result, [0, 35.728402, 104.601328, 54.777191])
+    assert_allclose(result, [0, 35.728402, 104.601328, 54.777191], atol=1e-4)
 
 
 def test_klucher_series_float():
@@ -120,29 +120,29 @@ def test_klucher_series():
     result = irradiance.klucher(40, 180, irrad_data['dhi'], irrad_data['ghi'],
                        ephem_data['apparent_zenith'],
                        ephem_data['azimuth'])
-    assert_allclose(result, [0, 37.446276, 109.209347, 56.965916])
+    assert_allclose(result, [0, 37.446276, 109.209347, 56.965916], atol=1e-4)
 
 
 def test_haydavies():
     result = irradiance.haydavies(40, 180, irrad_data['dhi'], irrad_data['dni'],
                          dni_et,
                          ephem_data['apparent_zenith'],
                          ephem_data['azimuth'])
-    assert_allclose(result, [0, 14.967008, 102.994862, 33.190865])
+    assert_allclose(result, [0, 14.967008, 102.994862, 33.190865], atol=1e-4)
 
 
 def test_reindl():
     result = irradiance.reindl(40, 180, irrad_data['dhi'], irrad_data['dni'],
                       irrad_data['ghi'], dni_et,
                       ephem_data['apparent_zenith'],
                       ephem_data['azimuth'])
-    assert_allclose(result, [np.nan, 15.730664, 104.131724, 34.166258])
+    assert_allclose(result, [np.nan, 15.730664, 104.131724, 34.166258], atol=1e-4)
 
 
 def test_king():
     result = irradiance.king(40, irrad_data['dhi'], irrad_data['ghi'],
                     ephem_data['apparent_zenith'])
-    assert_allclose(result, [0, 44.629352, 115.182626, 79.719855])
+    assert_allclose(result, [0, 44.629352, 115.182626, 79.719855], atol=1e-4)
 
 
 def test_perez():

pvlib/test/test_solarposition.py

@@ -36,6 +36,14 @@ def expected_solpos():
                          'apparent_elevation': 39.888378},
                         index=['2003-10-17T12:30:30Z'])
 
+@pytest.fixture()
+def expected_solpos_multi():
+    return pd.DataFrame({'elevation': [39.872046, 39.505196],
+                         'apparent_zenith': [50.111622, 50.478260],
+                         'azimuth': [194.340241, 194.311132],
+                         'apparent_elevation': [39.888378, 39.521740]},
+                        index=[['2003-10-17T12:30:30Z', '2003-10-18T12:30:30Z']])
+
 # the physical tests are run at the same time as the NREL SPA test.
 # pyephem reproduces the NREL result to 2 decimal places.
 # this doesn't mean that one code is better than the other.
@@ -63,7 +71,6 @@ def test_spa_c_physical_dst(expected_solpos):
     expected_solpos.index = times
     assert_frame_equal(expected_solpos, ephem_data[expected_solpos.columns])
 
-
 def test_spa_python_numpy_physical(expected_solpos):
     times = pd.date_range(datetime.datetime(2003,10,17,12,30,30),
                           periods=1, freq='D', tz=golden_mst.tz)
@@ -76,7 +83,6 @@ def test_spa_python_numpy_physical(expected_solpos):
     expected_solpos.index = times
     assert_frame_equal(expected_solpos, ephem_data[expected_solpos.columns])
 
-
 def test_spa_python_numpy_physical_dst(expected_solpos):
     times = pd.date_range(datetime.datetime(2003,10,17,13,30,30),
                           periods=1, freq='D', tz=golden.tz)
@@ -312,6 +318,30 @@ def test_get_solarposition_altitude(altitude, expected):
     assert_frame_equal(this_expected, ephem_data[this_expected.columns])
 
 
+@pytest.mark.parametrize(
+    "delta_t, method, expected", [
+    (None, 'nrel_numpy', expected_solpos_multi()),
+    (67.0, 'nrel_numpy', expected_solpos_multi()),
+    pytest.mark.xfail(raises=ValueError, reason = 'spa.calculate_deltat not implemented for numba yet')
+    ((None, 'nrel_numba', expected_solpos_multi())),
+    (67.0, 'nrel_numba', expected_solpos_multi())
+    ])
+def test_get_solarposition_deltat(delta_t, method, expected):
+    times = pd.date_range(datetime.datetime(2003,10,17,13,30,30),
+                          periods=2, freq='D', tz=golden.tz)
+    ephem_data = solarposition.get_solarposition(times, golden.latitude,
+                                                 golden.longitude,
+                                                 pressure=82000,
+                                                 delta_t=delta_t,
+                                                 temperature=11,
+                                                 method=method)
+    this_expected = expected.copy()
+    this_expected.index = times
+    this_expected = np.round(this_expected, 5)
+    ephem_data = np.round(ephem_data, 5)
+    assert_frame_equal(this_expected, ephem_data[this_expected.columns])
+
+
 def test_get_solarposition_no_kwargs(expected_solpos):
     times = pd.date_range(datetime.datetime(2003,10,17,13,30,30),
                           periods=1, freq='D', tz=golden.tz)