Merge branch 'main' into pr1725-continuation-shaded-fraction-mikofski

Author: echedey-ls

docs/sphinx/source/whatsnew/v0.10.4.rst

@@ -17,6 +17,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
@@ -37,8 +48,11 @@ 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`)
 * Mark A. Mikofski (:ghuser:`mikofski`)
-* Echedey Luis (:ghuser:`echedey-ls`)

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,21 +1387,23 @@ 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 <T>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.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)
 
 

pvlib/spa_c_files/SPA_NOTICE.md

@@ -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)

pvlib/tests/conftest.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