diff --git a/.stickler.yml b/.stickler.yml index 41aa301c06..dc360a6a7b 100644 --- a/.stickler.yml +++ b/.stickler.yml @@ -2,7 +2,7 @@ linters: flake8: python: 3 max-line-length: 79 - ignore: E201,E241,E226 + ignore: E201,E241,E226,W503,W504 files: ignore: - 'pvlib/_version.py' diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst index 45ef9ec50a..22e84d10a0 100644 --- a/docs/sphinx/source/api.rst +++ b/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 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 diff --git a/docs/sphinx/source/whatsnew/v0.7.0.rst b/docs/sphinx/source/whatsnew/v0.7.0.rst index cf2adfc9e6..a5c6f52e48 100644 --- a/docs/sphinx/source/whatsnew/v0.7.0.rst +++ b/docs/sphinx/source/whatsnew/v0.7.0.rst @@ -8,8 +8,11 @@ recommend all users of v0.6.3 upgrade to this release. **Python 2.7 support ended on June 1, 2019**. (:issue:`501`) +- Add a faster way to calculate the solar position (`spencer`) + Contributors ~~~~~~~~~~~~ * Mark Campanellli (:ghuser:`markcampanelli`) * Will Holmgren (:ghuser:`wholmgren`) +* Daniel Lassahn (:ghuser:`meteoDaniel`) \ No newline at end of file diff --git a/pvlib/solarposition.py b/pvlib/solarposition.py index d55737040b..ef2ca06eaa 100644 --- a/pvlib/solarposition.py +++ b/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,17 +19,16 @@ from imp import reload except ImportError: pass - import numpy as np import pandas as pd import warnings from pvlib import atmosphere -from pvlib.tools import datetime_to_djd, djd_to_datetime +from pvlib.tools import datetime_to_djd, djd_to_datetime, datetime_to_julian from pvlib._deprecation import deprecated - NS_PER_HR = 1.e9 * 3600. # nanoseconds per hour +JULIAN_YEARS = 365.2425 def get_solarposition(time, latitude, longitude, @@ -67,6 +67,9 @@ def get_solarposition(time, latitude, longitude, 'nrel_c' uses the NREL SPA C code [3]: :py:func:`spa_c` + 'spencer' uses the Spencer formula [4] :py:func:`spencer` + + temperature : float, default 12 Degrees C. @@ -81,6 +84,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 +118,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': + ephem_df = spencer(time, latitude, longitude) + else: raise ValueError('Invalid solar position method') @@ -495,7 +503,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 +1453,80 @@ 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(times, latitude, longitude): + """ + Calculate the solar position using a formulation by Spencer 1971/1972. + + Parameters + ---------- + times : 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), + 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, julian_2000 = datetime_to_julian(times) + julians_2000 = np.asarray(julians, dtype=np.float) - julian_2000 + + latitude_radians = np.radians(latitude) + day_time = (julians_2000 % 1) * 24 + + declination = np.array(declination_spencer71(times.dayofyear)) + + # Equation of time (difference between standard time and solar time). + eot = np.array(equation_of_time_spencer71(times.dayofyear)) + + # 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(latitude_radians) + + np.cos(declination) * np.cos(latitude_radians) * 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(latitude_radians) * np.sin(elevation) - np.sin(declination)) \ + / (np.cos(latitude_radians) * 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), + 'declination': declination, + 'equation_of_time': eot}, + index=times) + return result diff --git a/pvlib/spa.py b/pvlib/spa.py index 3a44a2efed..93e111e34d 100644 --- a/pvlib/spa.py +++ b/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. diff --git a/pvlib/test/test_solarposition.py b/pvlib/test/test_solarposition.py index 4a8c142dff..9052071573 100644 --- a/pvlib/test/test_solarposition.py +++ b/pvlib/test/test_solarposition.py @@ -754,8 +754,28 @@ def test_sun_rise_set_transit_geometric(expected_rise_set_spa, golden_mst): atol=np.abs(expected_transit_error).max()) -# put numba tests at end of file to minimize reloading +def test_get_solar_position_spencer(golden_mst): + """ test for the calculation based on spencer 1972 """ + times = pd.date_range(datetime.datetime(2003, 10, 17, 12, 30, 30), + periods=1, freq='D') + spencer_data = solarposition.spencer(times, + golden_mst.latitude, + golden_mst.longitude) + + expected_solpos = golden_mst.get_solarposition(times) + np.testing.assert_array_almost_equal(spencer_data.zenith.values, + expected_solpos.zenith.values, + decimal=0) + np.testing.assert_array_almost_equal(spencer_data.azimuth.values, + expected_solpos.azimuth.values, + decimal=0) + np.testing.assert_array_almost_equal(spencer_data.elevation.values, + expected_solpos.elevation.values, + decimal=0) + + +# put numba tests at end of file to minimize reloading @requires_numba def test_spa_python_numba_physical(expected_solpos, golden_mst): times = pd.date_range(datetime.datetime(2003, 10, 17, 12, 30, 30), diff --git a/pvlib/test/test_tools.py b/pvlib/test/test_tools.py index 42eaedca58..4a8fe73eb8 100644 --- a/pvlib/test/test_tools.py +++ b/pvlib/test/test_tools.py @@ -1,7 +1,14 @@ +import datetime + +import numpy as np +import pandas as pd import pytest from pvlib import tools +times = pd.date_range(start=datetime.datetime(2014, 6, 24), + end=datetime.datetime(2014, 6, 26), freq='15Min') + @pytest.mark.parametrize('keys, input_dict, expected', [ (['a', 'b'], {'a': 1, 'b': 2, 'c': 3}, {'a': 1, 'b': 2}), @@ -12,3 +19,21 @@ def test_build_kwargs(keys, input_dict, expected): kwargs = tools._build_kwargs(keys, input_dict) assert kwargs == expected + + +def test_datetime_to_julian(): + """ test transformation from datetime to julians """ + julians = tools.datetime_to_julian(pd.to_datetime(times)) + np.testing.assert_array_almost_equal(np.array(julians[:10]), + np.array([ + 2456832.5, + 2456832.5104166665, + 2456832.5208333335, + 2456832.53125, + 2456832.5416666665, + 2456832.5520833335, + 2456832.5625, + 2456832.5729166665, + 2456832.5833333335, + 2456832.59375]) + ) diff --git a/pvlib/tools.py b/pvlib/tools.py index cbd59a8541..065965a044 100644 --- a/pvlib/tools.py +++ b/pvlib/tools.py @@ -9,6 +9,10 @@ import pandas as pd import pytz +JULIAN_2000 = 2451545 +DAY_SECONDS = 60 * 60 * 24 +DT_2000 = pd.to_datetime(dt.datetime(2000, 1, 1)).tz_localize('UTC') + def cosd(angle): """ @@ -425,3 +429,30 @@ def _golden_sect_DataFrame(params, VL, VH, func): raise Exception("EXCEPTION:iterations exceeded maximum (50)") return func(df, 'V1'), df['V1'] + + +def datetime_to_julian(times): + """ + Transforms pd.DateTimeIndex to Julian days + + Parameters + ---------- + times : pandas.DatetimeIndex + Corresponding timestamps, must be localized to the timezone for the + ``latitude`` and ``longitude`` + Returns + ------- + Float64Index + The float index contains julian dates + Float + julian 2000 in UTC referenced to local time + """ + dt_2000 = DT_2000.tz_convert(times.tzinfo) + hours_difference = (DT_2000.tz_localize(None) - + dt_2000.tz_localize(None)).seconds / 3600 + julian_2000 = JULIAN_2000 - (hours_difference/24) + delta = times - dt_2000 + delta_julians = (delta.seconds + delta.microseconds / 1e6) + return ( + julian_2000 + delta.days + delta_julians / DAY_SECONDS + ), julian_2000