Skip to content

Fix various solar position calculations when DatetimeIndex.unit is not ns #1948

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 5 commits into from
Feb 23, 2024
Merged

Fix various solar position calculations when DatetimeIndex.unit is not ns #1948

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
db9b9d2
initial bugfixes and tests
kandersolar Jan 14, 2024
a7cc3b4
whatsnew
kandersolar Jan 14, 2024
356f5b7
compatibility with older pandas
kandersolar Jan 17, 2024
6f09d81
make tests make sense for pandas < and > 2
kandersolar Jan 17, 2024
2ea61bc
Merge branch 'main' into solpos-microseconds-fix
kandersolar Feb 23, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 15 additions & 0 deletions docs/sphinx/source/whatsnew/v0.10.4.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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`)
45 changes: 27 additions & 18 deletions pvlib/solarposition.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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',
Expand Down Expand Up @@ -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):
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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')
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should we warn about this localization? Maybe just a note. The function promises hour angle in the time zone for the longitude.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In case it's not obvious, this localization retains the existing behavior of assuming naive inputs to be in UTC; the rewrite for the unit issue just happens to make it more explicit. See for v0.9.5:

In [1]: import pvlib  # v0.9.5
   ...: import pandas as pd
   ...: for tz in [None, 'utc', 'Etc/GMT+5']:
   ...:     times = pd.date_range('2019-06-01 12:00', freq='h', periods=2, tz=tz)
   ...:     eot = pvlib.solarposition.equation_of_time_spencer71(times.dayofyear)
   ...:     ha = pvlib.solarposition.hour_angle(times, -80, eot)
   ...:     print(f'hour angles, {tz=}:', ha)
   ...:
hour angles, tz=None: [-79.36192279 -64.36192279]
hour angles, tz='utc': [-79.36192279 -64.36192279]
hour angles, tz='Etc/GMT+5': [-4.36192279 10.63807721]

The docstring says that times must be localized (Corresponding timestamps, must be localized to the timezone for the longitude.), so I guess naive times are technically not supported. Other solar position functions assume UTC for naive inputs, and don't warn, so I think I'm inclined to keep that behavior here as well.

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)

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


Expand Down
5 changes: 5 additions & 0 deletions pvlib/tests/conftest.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
118 changes: 116 additions & 2 deletions pvlib/tests/test_solarposition.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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),
Expand Down Expand Up @@ -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
Expand Down