diff --git a/docs/sphinx/source/whatsnew/v0.10.4.rst b/docs/sphinx/source/whatsnew/v0.10.4.rst index 5c529f6383..997ca4c249 100644 --- a/docs/sphinx/source/whatsnew/v0.10.4.rst +++ b/docs/sphinx/source/whatsnew/v0.10.4.rst @@ -12,6 +12,17 @@ Enhancements Bug fixes ~~~~~~~~~ +* Fixed an error in solar position calculations when using + :py:class:`pandas.DatetimeIndex` with ``unit`` other than ``'ns'`` (:issue:`1932`). + The following functions were affected: + + - :py:class:`~pvlib.modelchain.ModelChain` and :py:func:`~pvlib.solarposition.get_solarposition` with the ``nrel_numpy`` and ``nrel_numba`` methods + - :py:func:`~pvlib.solarposition.spa_python` + - :py:func:`~pvlib.solarposition.sun_rise_set_transit_spa` + - :py:func:`~pvlib.solarposition.nrel_earthsun_distance` + - :py:func:`~pvlib.solarposition.hour_angle` + - :py:func:`~pvlib.solarposition.sun_rise_set_transit_geometric` + * :py:class:`~pvlib.modelchain.ModelChain` now raises a more useful error when ``temperature_model_parameters`` are specified on the passed ``system`` instead of on its ``arrays``. (:issue:`1759`). * :py:func:`pvlib.irradiance.ghi_from_poa_driesse_2023` now correctly makes use @@ -32,6 +43,10 @@ Requirements Contributors ~~~~~~~~~~~~ +* Patrick Sheehan (:ghuser:`patricksheehan`) +* Echedey Luis (:ghuser:`echedey-ls`) +* Kevin Anderson (:ghuser:`kandersolar`) * Cliff Hansen (:ghuser:`cwhanse`) * :ghuser:`matsuobasho` * Adam R. Jensen (:ghuser:`AdamRJensen`) +* Peter Dudfield (:ghuser:`peterdudfield`) diff --git a/pvlib/solarposition.py b/pvlib/solarposition.py index 08023daf91..af59184727 100644 --- a/pvlib/solarposition.py +++ b/pvlib/solarposition.py @@ -28,9 +28,6 @@ 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, method='nrel_numpy', @@ -273,6 +270,19 @@ def _spa_python_import(how): return spa +def _datetime_to_unixtime(dtindex): + # convert a pandas datetime index to unixtime, making sure to handle + # different pandas units (ns, us, etc) and time zones correctly + if dtindex.tz is not None: + # epoch is 1970-01-01 00:00 UTC, but we need to match the input tz + # for compatibility with older pandas versions (e.g. v1.3.5) + epoch = pd.Timestamp("1970-01-01", tz="UTC").tz_convert(dtindex.tz) + else: + epoch = pd.Timestamp("1970-01-01") + + return np.array((dtindex - epoch) / pd.Timedelta("1s")) + + def spa_python(time, latitude, longitude, altitude=0, pressure=101325, temperature=12, delta_t=67.0, atmos_refract=None, how='numpy', numthreads=4): @@ -365,7 +375,7 @@ def spa_python(time, latitude, longitude, except (TypeError, ValueError): time = pd.DatetimeIndex([time, ]) - unixtime = np.array(time.view(np.int64)/10**9) + unixtime = _datetime_to_unixtime(time) spa = _spa_python_import(how) @@ -444,7 +454,7 @@ def sun_rise_set_transit_spa(times, latitude, longitude, how='numpy', # must convert to midnight UTC on day of interest utcday = pd.DatetimeIndex(times.date).tz_localize('UTC') - unixtime = np.array(utcday.view(np.int64)/10**9) + unixtime = _datetime_to_unixtime(utcday) spa = _spa_python_import(how) @@ -1000,7 +1010,7 @@ def nrel_earthsun_distance(time, how='numpy', delta_t=67.0, numthreads=4): except (TypeError, ValueError): time = pd.DatetimeIndex([time, ]) - unixtime = np.array(time.view(np.int64)/10**9) + unixtime = _datetime_to_unixtime(time) spa = _spa_python_import(how) @@ -1377,11 +1387,13 @@ def hour_angle(times, longitude, equation_of_time): equation_of_time_spencer71 equation_of_time_pvcdrom """ - naive_times = times.tz_localize(None) # naive but still localized # hours - timezone = (times - normalized_times) - (naive_times - times) - hrs_minus_tzs = 1 / NS_PER_HR * ( - 2 * times.view(np.int64) - times.normalize().view(np.int64) - - naive_times.view(np.int64)) + if times.tz is None: + times = times.tz_localize('utc') + tzs = np.array([ts.utcoffset().total_seconds() for ts in times]) / 3600 + + hrs_minus_tzs = (times - times.normalize()) / pd.Timedelta('1h') - tzs + # ensure array return instead of a version-dependent pandas Index return np.asarray( 15. * (hrs_minus_tzs - 12.) + longitude + equation_of_time / 4.) @@ -1389,9 +1401,9 @@ def hour_angle(times, longitude, equation_of_time): 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.view(np.int64) - times.view(np.int64)) + if times.tz is None: + times = times.tz_localize('utc') + tzs = np.array([ts.utcoffset().total_seconds() for ts in times]) / 3600 hours = (hourangle - longitude - equation_of_time / 4.) / 15. + 12. + tzs return np.asarray(hours) @@ -1405,16 +1417,13 @@ def _local_times_from_hours_since_midnight(times, hours): # normalize local, naive times to previous midnight and add the hours until # sunrise, sunset, and transit return pd.DatetimeIndex( - (naive_times.normalize().view(np.int64) + - (hours * NS_PER_HR).astype(np.int64)).astype('datetime64[ns]'), - tz=tz_info) + naive_times.normalize() + pd.to_timedelta(hours, unit='h'), 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.view(np.int64) - times.normalize().view(np.int64)) + hrs = (times - times.normalize()) / pd.Timedelta('1h') return np.array(hrs) diff --git a/pvlib/tests/conftest.py b/pvlib/tests/conftest.py index d74dfb35b4..355d7f0a5a 100644 --- a/pvlib/tests/conftest.py +++ b/pvlib/tests/conftest.py @@ -174,6 +174,11 @@ def has_numba(): requires_pysam = pytest.mark.skipif(not has_pysam, reason="requires PySAM") +has_pandas_2_0 = Version(pd.__version__) >= Version("2.0.0") +requires_pandas_2_0 = pytest.mark.skipif(not has_pandas_2_0, + reason="requires pandas>=2.0.0") + + @pytest.fixture() def golden(): return Location(39.742476, -105.1786, 'America/Denver', 1830.14) diff --git a/pvlib/tests/test_solarposition.py b/pvlib/tests/test_solarposition.py index 17870de27e..63e32a4c3c 100644 --- a/pvlib/tests/test_solarposition.py +++ b/pvlib/tests/test_solarposition.py @@ -12,8 +12,9 @@ from pvlib.location import Location from pvlib import solarposition, spa -from .conftest import requires_ephem, requires_spa_c, requires_numba - +from .conftest import ( + requires_ephem, requires_spa_c, requires_numba, requires_pandas_2_0 +) # setup times and locations to be tested. times = pd.date_range(start=datetime.datetime(2014, 6, 24), @@ -717,6 +718,119 @@ def test_sun_rise_set_transit_geometric(expected_rise_set_spa, golden_mst): atol=np.abs(expected_transit_error).max()) +@pytest.mark.parametrize('tz', [None, 'utc', 'US/Eastern']) +def test__datetime_to_unixtime(tz): + # for pandas < 2.0 where "unit" doesn't exist in pd.date_range. note that + # unit of ns is the only option in pandas<2, and the default in pandas 2.x + times = pd.date_range(start='2019-01-01', freq='h', periods=3, tz=tz) + expected = times.view(np.int64)/10**9 + actual = solarposition._datetime_to_unixtime(times) + np.testing.assert_equal(expected, actual) + + +@requires_pandas_2_0 +@pytest.mark.parametrize('unit', ['ns', 'us', 's']) +@pytest.mark.parametrize('tz', [None, 'utc', 'US/Eastern']) +def test__datetime_to_unixtime_units(unit, tz): + kwargs = dict(start='2019-01-01', freq='h', periods=3) + times = pd.date_range(**kwargs, unit='ns', tz='UTC') + expected = times.view(np.int64)/10**9 + + times = pd.date_range(**kwargs, unit=unit, tz='UTC').tz_convert(tz) + actual = solarposition._datetime_to_unixtime(times) + np.testing.assert_equal(expected, actual) + + +@requires_pandas_2_0 +@pytest.mark.parametrize('method', [ + 'nrel_numpy', + 'ephemeris', + pytest.param('pyephem', marks=requires_ephem), + pytest.param('nrel_numba', marks=requires_numba), + pytest.param('nrel_c', marks=requires_spa_c), +]) +@pytest.mark.parametrize('tz', [None, 'utc', 'US/Eastern']) +def test_get_solarposition_microsecond_index(method, tz): + # https://github.com/pvlib/pvlib-python/issues/1932 + + kwargs = dict(start='2019-01-01', freq='H', periods=24, tz=tz) + + index_ns = pd.date_range(unit='ns', **kwargs) + index_us = pd.date_range(unit='us', **kwargs) + + sp_ns = solarposition.get_solarposition(index_ns, 40, -80, method=method) + sp_us = solarposition.get_solarposition(index_us, 40, -80, method=method) + + assert_frame_equal(sp_ns, sp_us, check_index_type=False) + + +@requires_pandas_2_0 +@pytest.mark.parametrize('tz', [None, 'utc', 'US/Eastern']) +def test_nrel_earthsun_distance_microsecond_index(tz): + # https://github.com/pvlib/pvlib-python/issues/1932 + + kwargs = dict(start='2019-01-01', freq='H', periods=24, tz=tz) + + index_ns = pd.date_range(unit='ns', **kwargs) + index_us = pd.date_range(unit='us', **kwargs) + + esd_ns = solarposition.nrel_earthsun_distance(index_ns) + esd_us = solarposition.nrel_earthsun_distance(index_us) + + assert_series_equal(esd_ns, esd_us, check_index_type=False) + + +@requires_pandas_2_0 +@pytest.mark.parametrize('tz', [None, 'utc', 'US/Eastern']) +def test_hour_angle_microsecond_index(tz): + # https://github.com/pvlib/pvlib-python/issues/1932 + + kwargs = dict(start='2019-01-01', freq='H', periods=24, tz=tz) + + index_ns = pd.date_range(unit='ns', **kwargs) + index_us = pd.date_range(unit='us', **kwargs) + + ha_ns = solarposition.hour_angle(index_ns, -80, 0) + ha_us = solarposition.hour_angle(index_us, -80, 0) + + np.testing.assert_equal(ha_ns, ha_us) + + +@requires_pandas_2_0 +@pytest.mark.parametrize('tz', ['utc', 'US/Eastern']) +def test_rise_set_transit_spa_microsecond_index(tz): + # https://github.com/pvlib/pvlib-python/issues/1932 + + kwargs = dict(start='2019-01-01', freq='H', periods=24, tz=tz) + + index_ns = pd.date_range(unit='ns', **kwargs) + index_us = pd.date_range(unit='us', **kwargs) + + rst_ns = solarposition.sun_rise_set_transit_spa(index_ns, 40, -80) + rst_us = solarposition.sun_rise_set_transit_spa(index_us, 40, -80) + + assert_frame_equal(rst_ns, rst_us, check_index_type=False) + + +@requires_pandas_2_0 +@pytest.mark.parametrize('tz', [None, 'utc', 'US/Eastern']) +def test_rise_set_transit_geometric_microsecond_index(tz): + # https://github.com/pvlib/pvlib-python/issues/1932 + + kwargs = dict(start='2019-01-01', freq='H', periods=24, tz=tz) + + index_ns = pd.date_range(unit='ns', **kwargs) + index_us = pd.date_range(unit='us', **kwargs) + + args = (40, -80, 0, 0) + rst_ns = solarposition.sun_rise_set_transit_geometric(index_ns, *args) + rst_us = solarposition.sun_rise_set_transit_geometric(index_us, *args) + + for times_ns, times_us in zip(rst_ns, rst_us): + # can't use a fancy assert function here since the units are different + assert all(times_ns == times_us) + + # put numba tests at end of file to minimize reloading @requires_numba