diff --git a/docs/sphinx/source/api.rst b/docs/sphinx/source/api.rst index 2a232bc3fc..2c99b1608d 100644 --- a/docs/sphinx/source/api.rst +++ b/docs/sphinx/source/api.rst @@ -79,6 +79,7 @@ calculations. solarposition.equation_of_time_spencer71 solarposition.equation_of_time_pvcdrom solarposition.hour_angle + solarposition.sunrise_sunset_transit_geometric Clear sky diff --git a/docs/sphinx/source/whatsnew/v0.6.1.rst b/docs/sphinx/source/whatsnew/v0.6.1.rst index e9cc42e51c..7633e2dda1 100644 --- a/docs/sphinx/source/whatsnew/v0.6.1.rst +++ b/docs/sphinx/source/whatsnew/v0.6.1.rst @@ -3,7 +3,7 @@ v0.6.1 (EST October, 2018) -------------------------- -This is a minor release. We recommend all users of 0.6.0 upgrade to this +This is a minor release. We recommend all users of v0.6.0 upgrade to this release. **Python 2.7 support will end on June 1, 2019**. Releases made after this @@ -16,23 +16,30 @@ API Changes * Deprecated ``tmy``, ``tmy.readtmy2`` and ``tmy.readtmy3``; they will be removed in v0.7. Use the new :py:func:`pvlib.iotools.read_tmy2` and :py:func:`pvlib.iotools.read_tmy3` instead. (:issue:`261`) -* Added kwarg `horizon` to :func:`~pvlib.solarposition.pyephem` and :func:`~pvlib.solarposition.calc_time` with default value `'+0:00'` +* Added keyword argument ``horizon`` to :func:`~pvlib.solarposition.pyephem` + and :func:`~pvlib.solarposition.calc_time` with default value ``'+0:00'`` Enhancements ~~~~~~~~~~~~ -* :func:`~pvlib.solarposition.rise_set_transit_ephem` returns sunrise, sunset and transit times using pyephem -* Created :py:func:`pvlib.iotools.read_srml` and :py:func:`pvlib.iotools.read_srml_month_from_solardat` - to read University of Oregon Solar Radiation Monitoring Laboratory data. (:issue:`589`) - +* :func:`~pvlib.solarposition.rise_set_transit_ephem` returns sunrise, sunset + and transit times using pyephem (:issue:`114`) +* Add geometric functions for sunrise, sunset, and sun transit times, + :func:`~pvlib.solarposition.sunrise_sunset_transit_geometric` (:issue:`114`) +* Created :py:func:`pvlib.iotools.read_srml` and + :py:func:`pvlib.iotools.read_srml_month_from_solardat` to read University of + Oregon Solar Radiation Monitoring Laboratory data. (:issue:`589`) + Bug fixes ~~~~~~~~~ * Fix when building documentation using Matplotlib 3.0 or greater. +* Fix and improve :func:`~pvlib.solarposition.hour_angle` (:issue:`598`) Testing ~~~~~~~ +* Add test for :func:`~pvlib.solarposition.hour_angle` (:issue:`597`) Contributors @@ -40,3 +47,4 @@ Contributors * Will Holmgren (:ghuser:`wholmgren`) * Cliff Hansen (:ghuser:`cwhanse`) * Leland Boeman (:ghuser:`lboeman`) +* Mark Mikofski (:ghuser:`mikofski`) diff --git a/pvlib/solarposition.py b/pvlib/solarposition.py index 7d3db9d949..06a3b3e617 100644 --- a/pvlib/solarposition.py +++ b/pvlib/solarposition.py @@ -25,6 +25,8 @@ from pvlib import atmosphere from pvlib.tools import datetime_to_djd, djd_to_datetime +NS_PER_HR = 1.e9 * 3600. # nanoseconds per hour + def get_solarposition(time, latitude, longitude, altitude=None, pressure=None, @@ -1175,7 +1177,7 @@ def declination_cooper69(dayofyear): return dec -def solar_azimuth_analytical(latitude, hour_angle, declination, zenith): +def solar_azimuth_analytical(latitude, hourangle, declination, zenith): """ Analytical expression of solar azimuth angle based on spherical trigonometry. @@ -1184,7 +1186,7 @@ def solar_azimuth_analytical(latitude, hour_angle, declination, zenith): ---------- latitude : numeric Latitude of location in radians. - hour_angle : numeric + hourangle : numeric Hour angle in the local solar time in radians. declination : numeric Declination of the sun in radians. @@ -1240,12 +1242,12 @@ def solar_azimuth_analytical(latitude, hour_angle, declination, zenith): # when NaN values occur in input, ignore and pass to output with np.errstate(invalid='ignore'): - sign_ha = np.sign(hour_angle) + sign_ha = np.sign(hourangle) - return (sign_ha * np.arccos(cos_azi) + np.pi) + return sign_ha * np.arccos(cos_azi) + np.pi -def solar_zenith_analytical(latitude, hour_angle, declination): +def solar_zenith_analytical(latitude, hourangle, declination): """ Analytical expression of solar zenith angle based on spherical trigonometry. @@ -1257,7 +1259,7 @@ def solar_zenith_analytical(latitude, hour_angle, declination): ---------- latitude : numeric Latitude of location in radians. - hour_angle : numeric + hourangle : numeric Hour angle in the local solar time in radians. declination : numeric Declination of the sun in radians. @@ -1288,7 +1290,7 @@ def solar_zenith_analytical(latitude, hour_angle, declination): declination_spencer71 declination_cooper69 hour_angle """ return np.arccos( - np.cos(declination) * np.cos(latitude) * np.cos(hour_angle) + + np.cos(declination) * np.cos(latitude) * np.cos(hourangle) + np.sin(declination) * np.sin(latitude) ) @@ -1300,7 +1302,8 @@ def hour_angle(times, longitude, equation_of_time): Parameters ---------- times : :class:`pandas.DatetimeIndex` - Corresponding timestamps, must be timezone aware. + Corresponding timestamps, must be localized to the timezone for the + ``longitude``. longitude : numeric Longitude in degrees equation_of_time : numeric @@ -1329,9 +1332,99 @@ def hour_angle(times, longitude, equation_of_time): """ naive_times = times.tz_localize(None) # naive but still localized # hours - timezone = (times - normalized_times) - (naive_times - times) - hrs_minus_tzs = 1 / (3600. * 1.e9) * ( + hrs_minus_tzs = 1 / NS_PER_HR * ( 2 * times.astype(np.int64) - times.normalize().astype(np.int64) - naive_times.astype(np.int64)) # ensure array return instead of a version-dependent pandas Index return np.asarray( 15. * (hrs_minus_tzs - 12.) + longitude + equation_of_time / 4.) + + +def _hour_angle_to_hours(times, hourangle, longitude, equation_of_time): + """converts hour angles in degrees to hours as a numpy array""" + naive_times = times.tz_localize(None) # naive but still localized + tzs = 1 / NS_PER_HR * ( + naive_times.astype(np.int64) - times.astype(np.int64)) + hours = (hourangle - longitude - equation_of_time / 4.) / 15. + 12. + tzs + return np.asarray(hours) + + +def _local_times_from_hours_since_midnight(times, hours): + """ + converts hours since midnight from an array of floats to localized times + """ + tz_info = times.tz # pytz timezone info + naive_times = times.tz_localize(None) # naive but still localized + # normalize local, naive times to previous midnight and add the hours until + # sunrise, sunset, and transit + return pd.DatetimeIndex( + (naive_times.normalize().astype(np.int64) + + (hours * NS_PER_HR).astype(np.int64)).astype('datetime64[ns]'), + tz=tz_info) + + +def _times_to_hours_after_local_midnight(times): + """convert local pandas datetime indices to array of hours as floats""" + times = times.tz_localize(None) + hrs = 1 / NS_PER_HR * ( + times.astype(np.int64) - times.normalize().astype(np.int64)) + return np.array(hrs) + + +def sunrise_sunset_transit_geometric(times, latitude, longitude, declination, + equation_of_time): + """ + Geometric calculation of solar sunrise, sunset, and transit. + + .. warning:: The geometric calculation assumes a circular earth orbit with + the sun as a point source at its center, and neglects the effect of + atmospheric refraction on zenith. The error depends on location and + time of year but is of order 10 minutes. + + 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 + declination : numeric + declination angle in radians at ``times`` + equation_of_time : numeric + difference in time between solar time and mean solar time in minutes + + Returns + ------- + sunrise : datetime + localized sunrise time + sunset : datetime + localized sunset time + transit : datetime + localized sun transit time + + References + ---------- + [1] J. A. Duffie and W. A. Beckman, "Solar Engineering of Thermal + Processes, 3rd Edition," J. Wiley and Sons, New York (2006) + + [2] Frank Vignola et al., "Solar And Infrared Radiation Measurements," + CRC Press (2012) + + """ + latitude_rad = np.radians(latitude) # radians + sunset_angle_rad = np.arccos(-np.tan(declination) * np.tan(latitude_rad)) + sunset_angle = np.degrees(sunset_angle_rad) # degrees + # solar noon is at hour angle zero + # so sunrise is just negative of sunset + sunrise_angle = -sunset_angle + sunrise_hour = _hour_angle_to_hours( + times, sunrise_angle, longitude, equation_of_time) + sunset_hour = _hour_angle_to_hours( + times, sunset_angle, longitude, equation_of_time) + transit_hour = _hour_angle_to_hours(times, 0, longitude, equation_of_time) + sunrise = _local_times_from_hours_since_midnight(times, sunrise_hour) + sunset = _local_times_from_hours_since_midnight(times, sunset_hour) + transit = _local_times_from_hours_since_midnight(times, transit_hour) + return sunrise, sunset, transit diff --git a/pvlib/test/test_solarposition.py b/pvlib/test/test_solarposition.py index efc424e8a1..c835a357a0 100644 --- a/pvlib/test/test_solarposition.py +++ b/pvlib/test/test_solarposition.py @@ -529,11 +529,13 @@ def test_get_solarposition_altitude(altitude, expected, golden): @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()) + (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, golden): times = pd.date_range(datetime.datetime(2003,10,17,13,30,30), @@ -690,8 +692,9 @@ def test_analytical_azimuth(): [ -30., -180., 5.], [ -30., 180., 10.]])) - zeniths = solarposition.solar_zenith_analytical(*test_angles.T) - azimuths = solarposition.solar_azimuth_analytical(*test_angles.T, zenith=zeniths) + zeniths = solarposition.solar_zenith_analytical(*test_angles.T) + azimuths = solarposition.solar_azimuth_analytical(*test_angles.T, + zenith=zeniths) assert not np.isnan(azimuths).any() @@ -718,3 +721,48 @@ def test_hour_angle(): # sunrise: 4 seconds, sunset: 48 seconds, transit: 0.2 seconds # but the differences may be due to other SPA input parameters assert np.allclose(hours, expected) + + +def test_sunrise_sunset_transit_geometric(expected_rise_set_spa, golden_mst): + """Test geometric calculations for sunrise, sunset, and transit times""" + times = expected_rise_set_spa.index + latitude = golden_mst.latitude + longitude = golden_mst.longitude + eot = solarposition.equation_of_time_spencer71(times.dayofyear) # minutes + decl = solarposition.declination_spencer71(times.dayofyear) # radians + sr, ss, st = solarposition.sunrise_sunset_transit_geometric( + times, latitude=latitude, longitude=longitude, declination=decl, + equation_of_time=eot) + # sunrise: 2015-01-02 07:26:39.763224487, 2015-08-02 05:04:35.688533801 + # sunset: 2015-01-02 16:41:29.951096777, 2015-08-02 19:09:46.597355085 + # transit: 2015-01-02 12:04:04.857160632, 2015-08-02 12:07:11.142944443 + test_sunrise = solarposition._times_to_hours_after_local_midnight(sr) + test_sunset = solarposition._times_to_hours_after_local_midnight(ss) + test_transit = solarposition._times_to_hours_after_local_midnight(st) + # convert expected SPA sunrise, sunset, transit to local datetime indices + expected_sunrise = pd.DatetimeIndex(expected_rise_set_spa.sunrise.values, + tz='UTC').tz_convert(golden_mst.tz) + expected_sunset = pd.DatetimeIndex(expected_rise_set_spa.sunset.values, + tz='UTC').tz_convert(golden_mst.tz) + expected_transit = pd.DatetimeIndex(expected_rise_set_spa.transit.values, + tz='UTC').tz_convert(golden_mst.tz) + # convert expected times to hours since midnight as arrays of floats + expected_sunrise = solarposition._times_to_hours_after_local_midnight( + expected_sunrise) + expected_sunset = solarposition._times_to_hours_after_local_midnight( + expected_sunset) + expected_transit = solarposition._times_to_hours_after_local_midnight( + expected_transit) + # geometric time has about 4-6 minute error compared to SPA sunset/sunrise + expected_sunrise_error = np.array( + [0.07910089555555544, 0.06908014805555496]) # 4.8[min], 4.2[min] + expected_sunset_error = np.array( + [-0.1036246955555562, -0.06983406805555603]) # -6.2[min], -4.2[min] + expected_transit_error = np.array( + [-0.011150788888889096, 0.0036508177777765383]) # -40[sec], 13.3[sec] + assert np.allclose(test_sunrise, expected_sunrise, + atol=np.abs(expected_sunrise_error).max()) + assert np.allclose(test_sunset, expected_sunset, + atol=np.abs(expected_sunset_error).max()) + assert np.allclose(test_transit, expected_transit, + atol=np.abs(expected_transit_error).max())