Solar position values are incorrect when using DatetimeIndex with unit='us'

[patricksheehan]

Describe the bug
Solar elevation values are wrong when using the default nrel_numpy method

To Reproduce
Steps to reproduce the behavior:

import pvlib

# Example usage:
latitude = 37.7749  # Replace with the desired latitude
longitude = -122.4194  # Replace with the desired longitude

# Create a sample DatetimeIndex (replace with your own datetime index)
date_index = pd.date_range(start='2023-01-01', end='2023-01-05', freq='H')

correct_results = pvlib.solarposition.get_solarposition(
    time=date_index,
    latitude=37.804363,
    longitude=-122.271111,
    method="pyephem"
)

incorrect_results = pvlib.solarposition.get_solarposition(
    time=date_index,
    latitude=37.804363,
    longitude=-122.271111,
) 

Expected behavior
solar altitude should be on a negative/positive cycle throughout a day

Screenshots

Versions:

• pvlib.__version__: 0.10.2
• pandas.__version__: 2.1.4
• python: 3.10.13

Additional context
Add any other context about the problem here.

There's obviously a workaround, but this was a surprising error, is this known?

[cwhanse]

I get the same results in both cases.

If the elevation isn't lined up with time, then the datetimes need to be localized with the correct timezone.

[patricksheehan]

Interesting, FWIW this is what the code looks like which generated the error on my end, and you can see the elevation data is off (no negatives throughout a day). Not sure it's worth figuring out reconstructing this exact index, but in case you see why this index might be different:

[kandersolar]

I think the problem is that the time index in that screenshot has dtype datetime64[us, US/Pacific], but spa_python is assuming the input is represented in nanoseconds:

pvlib-python/pvlib/solarposition.py

Line 369 in 8057393

unixtime = np.array(time.view(np.int64)/10**9)

Here is code to recreate the issue from scratch:

import pandas as pd
import pvlib

for unit in ['ns', 'us']:
    index = pd.date_range('2023-12-09', freq='H', periods=6, tz='US/Pacific', unit=unit)
    sp = pvlib.solarposition.get_solarposition(time=index, latitude=37.804363, longitude=-122.271111)
    print(f'Results for {unit=}:')
    print(sp['elevation'])

Results for unit='ns':
2023-12-09 00:00:00-08:00   -74.994002
2023-12-09 01:00:00-08:00   -70.408287
2023-12-09 02:00:00-08:00   -60.497812
2023-12-09 03:00:00-08:00   -49.067690
2023-12-09 04:00:00-08:00   -37.255747
2023-12-09 05:00:00-08:00   -25.482336
Freq: H, Name: elevation, dtype: float64
Results for unit='us':
2023-12-09 00:00:00-08:00    13.682861
2023-12-09 01:00:00-08:00    13.692004
2023-12-09 02:00:00-08:00    13.701146
2023-12-09 03:00:00-08:00    13.710287
2023-12-09 04:00:00-08:00    13.719425
2023-12-09 05:00:00-08:00    13.728563
Freq: H, Name: elevation, dtype: float64

@patricksheehan, just to be sure, when you run the original reproducer code in #1932 (comment), do you observe the issue? Or is it only with that us index in the notebook?

[cwhanse]

@kandersolar unit was introduced at pandas 2.0.0, I checked this with 1.5.3, maybe that's why I didn't reproduce the error.

[kandersolar]

Well, I think ns is still the default unit in pandas 2.x. The only way I know how to get an index with us units is to explicitly request it (either specifying unit='us' when creating the index or converting an existing index with .astype('datetime64[us]')). I'm curious how the index in @patricksheehan's example came to have us units :)

1e9 and 10**9 show up in several places in solarposition.py, so we might need to fix this problem in a few different places.

[patricksheehan]

@kandersolar thanks for looking at that, you're totally right that the us units is the culprit, the OG code I posted works for me as well (apologies for the bad test code, I always settle at 50% effort for reporting weird things vs not reporting at all).

I'm also honestly not sure how the us unit got in there, this is from a client's pickled DataFrame so who knows ;).

[cwhanse]

It's a good heads up for an issue that needs fixed, I don't think we want to limit pvlib to pandas < 2.0.0

[kandersolar]

Another instance of this bug in the wild (pandas 2.0.2): https://groups.google.com/g/pvlib-python/c/xDcgxTwmGTA

[echedey-ls]

I'd like to do my bit with this issue since I've also experienced this looking to optimize memory usage using rougher precision with datetime[s].

Indeed, on the Internet the // 10**9 solution is clearly favoured, and it really makes sense, since current pandas.DatetimeIndex docs state the following as kandersolar noted before.

dtype : numpy.dtype or DatetimeTZDtype or str, default None
    Note that the only NumPy dtype allowed is datetime64[ns].

1. This would be a reason to consider this a not-PVLIB issue. I'm not really in favour of that.

2. A possible solution that I haven't tried yet could be:

   1. Converting to datetime.datetime with pd.DatetimeIndex.to_pydatetime. It is not affected by the underlying dtype.
   2. Follow this StackOverflow post. Or this one, mentioned in the former.

3. Another possibility is to warn about the underlying type of pd.DatetimeIndex in the .get_solarposition docs.

4. Regarding doing the conversion manually in PVLIB, it is feasible. I've checked Pandas is capable of doing the right introspection of the dtype for non-ns data types:

>>> import pandas as pd
>>> from importlib.metadata import version
>>> pd.DatetimeIndex(['2023-09-06']).dtype
dtype('<M8[ns]')
>>> pd.DatetimeIndex(['2023-09-06'], dtype='datetime64[ns]').dtype
dtype('<M8[ns]')
>>> pd.DatetimeIndex(['2023-09-06'], dtype='datetime64[s]').dtype 
dtype('<M8[s]')
>>> pd.DatetimeIndex(pd.date_range('2023-09-03', '2023-09-06', unit='s')).dtype
dtype('<M8[s]')
>>> pd.DatetimeIndex(pd.date_range('2023-09-03', '2023-09-06', unit='s')).to_pydatetime()
array([datetime.datetime(2023, 9, 3, 0, 0),
       datetime.datetime(2023, 9, 4, 0, 0),
       datetime.datetime(2023, 9, 5, 0, 0),
       datetime.datetime(2023, 9, 6, 0, 0)], dtype=object)
>>> version("pandas")
'2.0.2'

I advocate for doing no. 2. .to_pydatetime exists since Pandas 1.0.0

I may work on this on the future, but if someone else wants, I encourage doing so.

[kandersolar]

Thanks @echedey-ls! I think we really need to have pvlib work with us and other units.

The pandas docs suggest this method of calculating unix time, which seems to be working well in my tests:

(stamps - pd.Timestamp("1970-01-01")) // pd.Timedelta("1s")

I have a PR in the works with this approach. The main thing is coming up with comprehensive tests...

[echedey-ls]

Oh, wow, I didn't find that earlier! That's incredibly Pythonic.

Good luck with the tests, that's the worst part.