diff --git a/docs/sphinx/source/whatsnew/v0.4.2.txt b/docs/sphinx/source/whatsnew/v0.4.2.txt index 9569a42db8..0b12a564ed 100644 --- a/docs/sphinx/source/whatsnew/v0.4.2.txt +++ b/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 diff --git a/pvlib/solarposition.py b/pvlib/solarposition.py index 2e124c33f0..276a469650 100644 --- a/pvlib/solarposition.py +++ b/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) diff --git a/pvlib/spa.py b/pvlib/spa.py index 80aae7c24b..770ad110e8 100644 --- a/pvlib/spa.py +++ b/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 diff --git a/pvlib/test/test_irradiance.py b/pvlib/test/test_irradiance.py index 573ae40f57..83750cb6e2 100644 --- a/pvlib/test/test_irradiance.py +++ b/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,7 +120,7 @@ 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(): @@ -128,7 +128,7 @@ def test_haydavies(): 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(): @@ -136,13 +136,13 @@ def test_reindl(): 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(): diff --git a/pvlib/test/test_solarposition.py b/pvlib/test/test_solarposition.py index af39e77332..c4dc5c1200 100644 --- a/pvlib/test/test_solarposition.py +++ b/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) diff --git a/pvlib/test/test_spa.py b/pvlib/test/test_spa.py index 1962bca5c0..740948312a 100644 --- a/pvlib/test/test_spa.py +++ b/pvlib/test/test_spa.py @@ -77,7 +77,19 @@ theta0 = 90 - e0 Gamma = 14.340241 Phi = 194.340241 - +year = 1985 +month = 2 +year_array = np.array([-499, 500, 1000, 1500, 1800, 1900, 1950, 1970, 1985, 1990, 2000, 2005]) +month_array = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]) +dt_actual = 54.413442486 +dt_actual_array = np.array([1.7184831e+04, 5.7088051e+03, 1.5730419e+03, + 1.9801820e+02, 1.3596506e+01, -2.1171894e+00, + 2.9289261e+01, 4.0824887e+01, 5.4724581e+01, + 5.7426651e+01, 6.4108015e+01, 6.5038015e+01]) +mix_year_array = np.full((10), year) +mix_month_array = np.full((10), month) +mix_year_actual = np.full((10), dt_actual) +mix_month_actual = mix_year_actual class SpaBase(object): """Test functions common to numpy and numba spa""" @@ -330,6 +342,19 @@ def test_earthsun_distance(self): result = self.spa.earthsun_distance(unixtimes, 64.0, 1) assert_almost_equal(R, result, 6) + def test_calculate_deltat(self): + result_mix_year = self.spa.calculate_deltat(mix_year_array, month) + assert_almost_equal(mix_year_actual, result_mix_year) + + result_mix_month = self.spa.calculate_deltat(year, mix_month_array) + assert_almost_equal(mix_month_actual, result_mix_month) + + result_array = self.spa.calculate_deltat(year_array, month_array) + assert_almost_equal(dt_actual_array, result_array, 3) + + result_scalar = self.spa.calculate_deltat(year,month) + assert_almost_equal(dt_actual, result_scalar) + class NumpySpaTest(unittest.TestCase, SpaBase): """Import spa without compiling to numba then run tests""" @classmethod