Perform dayofyear-based calculations according to UTC, not local time (#2055)

* Add clearsky.ipynb (Remove before merge)

* Fix clearsky.lookup_linke_turbidity to be timezone-aware

* Fix irradiance.get_extra_radiation to be timezone-aware

* Rerun and clear errors in clearsky.ipynb

* Remove clearsky.ipynb

* Update what's new

* Add docstring Notes to lookup_linke_turbidity

* Add tools._pandas_to_utc

* Update what's new

Co-authored-by: Cliff Hansen <cwhanse@sandia.gov>

* Add docstring Notes to _pandas_to_doy

* Revert _datetimelike_scalar_to_datetimeindex

* Revert tests.test_irradiance.times

* Replace lowercase utc with uppercase

* Fix clearsky lookup_linke_turbidity and _interpolate_turbidity

* Fix linter

* Fix clearsky._interpolate_turbidity to keep timezone

* Revert irradiance.py

* Clarify UTC input of ephem.Date

* Fix tz-naive delta_t of dni_extra nrel

* Fix tz-naive delta_t of get_solarposition nrel_numpy

* Avoid future midnight of location.get_sun_rise_set_transit spa

* Fix tz-naive delta_t of location.get_sun_rise_set_transit spa

* Fix local midnight of hour_angle and sun_rise_set_transit

* Fix dayofyear of solarposition decl and eot, i.e. clearsky.bird

* Fix dayofyear of spectrum.spectrl2 dni_extra spencer

* Fix monthly SeasonalTiltMount orientation

* Fix year and dayofyear of solrad

* Fix flake8

* resolve whatsnew

* don't change v0.11.0

* minor adjustment to test_clearsky

* undo edit to get_solrad

* undo edit to get_solrad

* Update pvlib/clearsky.py

* Update docs/sphinx/source/whatsnew/v0.11.1.rst

* Update docs/sphinx/source/whatsnew/v0.11.1.rst

* Update docs/examples/irradiance-transposition/plot_seasonal_tilt.py

Co-authored-by: Kevin Anderson <kevin.anderso@gmail.com>

* Revert benchmarks/benchmarks/solarposition.py

* Revert test_lookup_linke_turbidity_nointerp

* Revert solrad.py

* Revert plot_seasonal_tilt.py except notes

* Update test_sun_rise_set_transit_spa with delta_t

* Update test_nrel_earthsun_distance with delta_t

* Update test_sun_rise_set_transit_geometric with naive tz

* Update test_hour_angle with inversion and naive tz

* Remove _is_leap_year from _interpolate_turbidity

* Update docs/sphinx/source/whatsnew/v0.11.1.rst

---------

Co-authored-by: Cliff Hansen <cwhanse@sandia.gov>
Co-authored-by: Kevin Anderson <kevin.anderso@gmail.com>

Author: 3 people

docs/examples/irradiance-transposition/plot_seasonal_tilt.py

@@ -32,6 +32,8 @@ class SeasonalTiltMount(pvsystem.AbstractMount):
     surface_azimuth: float = 180.0
 
     def get_orientation(self, solar_zenith, solar_azimuth):
+        # note: determining tilt based on month may produce different
+        # results depending on the time zone of the timestamps
         tilts = [self.monthly_tilts[m-1] for m in solar_zenith.index.month]
         return pd.DataFrame({
             'surface_tilt': tilts,

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

@@ -39,7 +39,19 @@ Enhancements
 
 Bug fixes
 ~~~~~~~~~
-
+* To prevent simulation output from differing slightly based on the time zone
+  of the time stamps, models that use day of year for sun position and
+  irradiance calculations now determine the day of year according to the UTC
+  equivalent of the specified time stamps.  The following functions are affected:
+  :py:func:`~pvlib.clearsky.lookup_linke_turbidity`,
+  :py:func:`~pvlib.irradiance.get_extra_radiation`,
+  :py:func:`~pvlib.irradiance.disc`,
+  :py:func:`~pvlib.irradiance.dirint`,
+  :py:func:`~pvlib.spectrum.spectrl2`. (:issue:`2054`, :pull:`2055`)
+* :py:func:`~pvlib.solarposition.hour_angle` and
+  :py:func:`~pvlib.solarposition.sun_rise_set_transit_geometric` now raise
+  ``ValueError`` when given timezone-naive inputs, instead of assuming UTC.
+  (:pull:`2055`)
 
 Testing
 ~~~~~~~
@@ -70,6 +82,7 @@ Requirements
 Contributors
 ~~~~~~~~~~~~
 * Echedey Luis (:ghuser:`echedey-ls`)
+* Yunho Kee (:ghuser:`yhkee0404`)
 * Chris Deline (:ghuser:`cdeline`)
 * Ioannis Sifnaios (:ghuser:`IoannisSifnaios`)
 * Leonardo Micheli (:ghuser:`lmicheli`)

docs/tutorials/irradiance.ipynb

@@ -16,7 +16,8 @@
     "This tutorial requires pvlib >= 0.6.0.\n",
     "\n",
     "Authors:\n",
-    "* Will Holmgren (@wholmgren), University of Arizona. July 2014, April 2015, July 2015, March 2016, July 2016, February 2017, August 2018."
+    "* Will Holmgren (@wholmgren), University of Arizona. July 2014, April 2015, July 2015, March 2016, July 2016, February 2017, August 2018.\n",
+    "* Yunho Kee (@yhkee0404), Haezoom. Jun 2024."
    ]
   },
   {
@@ -393,6 +394,7 @@
    "source": [
     "tus = pvlib.location.Location(32.2, -111, 'US/Arizona', 700, 'Tucson')\n",
     "times = pd.date_range(start='2016-01-01', end='2016-01-02', freq='1min', tz=tus.tz)\n",
+    "times_utc = times.tz_convert('UTC')\n",
     "ephem_data = tus.get_solarposition(times)\n",
     "irrad_data = tus.get_clearsky(times)\n",
     "irrad_data.plot()\n",
@@ -810,7 +812,7 @@
     "                                            ephem_data['apparent_zenith'], ephem_data['azimuth'])\n",
     "klucher_diffuse.plot(label='klucher diffuse')\n",
     "\n",
-    "dni_et = pvlib.irradiance.get_extra_radiation(times.dayofyear)\n",
+    "dni_et = pvlib.irradiance.get_extra_radiation(times_utc.dayofyear)\n",
     "reindl_diffuse = pvlib.irradiance.reindl(surf_tilt, surf_az, \n",
     "                                          irrad_data['dhi'], irrad_data['dni'], irrad_data['ghi'], dni_et,\n",
     "                                          ephem_data['apparent_zenith'], ephem_data['azimuth'])\n",
@@ -858,7 +860,7 @@
     "                                            ephem_data['apparent_zenith'], ephem_data['azimuth'])\n",
     "klucher_diffuse.plot(label='klucher diffuse')\n",
     "\n",
-    "dni_et = pvlib.irradiance.get_extra_radiation(times.dayofyear)\n",
+    "dni_et = pvlib.irradiance.get_extra_radiation(times_utc.dayofyear)\n",
     "reindl_diffuse = pvlib.irradiance.reindl(surf_tilt, surf_az, \n",
     "                                          irrad_data['dhi'], irrad_data['dni'], irrad_data['ghi'], dni_et,\n",
     "                                          ephem_data['apparent_zenith'], ephem_data['azimuth'])\n",
@@ -913,7 +915,7 @@
     "                                            ephem_data['apparent_zenith'], ephem_data['azimuth'])\n",
     "klucher_diffuse.plot(label='klucher diffuse')\n",
     "\n",
-    "dni_et = pvlib.irradiance.get_extra_radiation(times.dayofyear)\n",
+    "dni_et = pvlib.irradiance.get_extra_radiation(times_utc.dayofyear)\n",
     "\n",
     "haydavies_diffuse = pvlib.irradiance.haydavies(surf_tilt, surf_az, \n",
     "                                                irrad_data['dhi'], irrad_data['dni'], dni_et,\n",
@@ -967,7 +969,7 @@
     "                                            ephem_data['apparent_zenith'], ephem_data['azimuth'])\n",
     "klucher_diffuse.plot(label='klucher diffuse')\n",
     "\n",
-    "dni_et = pvlib.irradiance.get_extra_radiation(times.dayofyear)\n",
+    "dni_et = pvlib.irradiance.get_extra_radiation(times_utc.dayofyear)\n",
     "\n",
     "haydavies_diffuse = pvlib.irradiance.haydavies(surf_tilt, surf_az, \n",
     "                                                irrad_data['dhi'], irrad_data['dni'], dni_et,\n",
@@ -1028,7 +1030,7 @@
     "                                            ephem_data['apparent_zenith'], ephem_data['azimuth'])\n",
     "klucher_diffuse.plot(label='klucher diffuse')\n",
     "\n",
-    "dni_et = pvlib.irradiance.get_extra_radiation(times.dayofyear)\n",
+    "dni_et = pvlib.irradiance.get_extra_radiation(times_utc.dayofyear)\n",
     "\n",
     "haydavies_diffuse = pvlib.irradiance.haydavies(surf_tilt, surf_az, \n",
     "                                                irrad_data['dhi'], irrad_data['dni'], dni_et,\n",
@@ -1067,7 +1069,7 @@
     "sun_az = ephem_data['azimuth']\n",
     "DNI = irrad_data['dni']\n",
     "DHI = irrad_data['dhi']\n",
-    "DNI_ET = pvlib.irradiance.get_extra_radiation(times.dayofyear)\n",
+    "DNI_ET = pvlib.irradiance.get_extra_radiation(times_utc.dayofyear)\n",
     "AM = pvlib.atmosphere.get_relative_airmass(sun_zen)\n",
     "\n",
     "surf_tilt = 32\n",
@@ -1316,7 +1318,7 @@
     "sun_az = ephem_data['azimuth']\n",
     "DNI = irrad_data['dni']\n",
     "DHI = irrad_data['dhi']\n",
-    "DNI_ET = pvlib.irradiance.get_extra_radiation(times.dayofyear)\n",
+    "DNI_ET = pvlib.irradiance.get_extra_radiation(times_utc.dayofyear)\n",
     "AM = pvlib.atmosphere.get_relative_airmass(sun_zen)\n",
     "\n",
     "surf_tilt = 32\n",
@@ -1330,7 +1332,7 @@
     "                                            ephem_data['apparent_zenith'], ephem_data['azimuth'])\n",
     "klucher_diffuse.plot(label='klucher diffuse')\n",
     "\n",
-    "dni_et = pvlib.irradiance.get_extra_radiation(times.dayofyear)\n",
+    "dni_et = pvlib.irradiance.get_extra_radiation(times_utc.dayofyear)\n",
     "\n",
     "haydavies_diffuse = pvlib.irradiance.haydavies(surf_tilt, surf_az, \n",
     "                                                irrad_data['dhi'], irrad_data['dni'], dni_et,\n",
@@ -1371,7 +1373,7 @@
     "sun_az = ephem_data['azimuth']\n",
     "DNI = irrad_data['dni']\n",
     "DHI = irrad_data['dhi']\n",
-    "DNI_ET = pvlib.irradiance.get_extra_radiation(times.dayofyear)\n",
+    "DNI_ET = pvlib.irradiance.get_extra_radiation(times_utc.dayofyear)\n",
     "AM = pvlib.atmosphere.get_relative_airmass(sun_zen)\n",
     "\n",
     "surf_tilt = 32\n",
@@ -1385,7 +1387,7 @@
     "                                            ephem_data['apparent_zenith'], ephem_data['azimuth'])\n",
     "klucher_diffuse.plot(label='klucher diffuse')\n",
     "\n",
-    "dni_et = pvlib.irradiance.get_extra_radiation(times.dayofyear)\n",
+    "dni_et = pvlib.irradiance.get_extra_radiation(times_utc.dayofyear)\n",
     "\n",
     "haydavies_diffuse = pvlib.irradiance.haydavies(surf_tilt, surf_az, \n",
     "                                                irrad_data['dhi'], irrad_data['dni'], dni_et,\n",

pvlib/clearsky.py

@@ -168,6 +168,13 @@ def lookup_linke_turbidity(time, latitude, longitude, filepath=None,
     Returns
     -------
     turbidity : Series
+
+    Notes
+    -----
+    Linke turbidity is obtained from a file of historical monthly averages.
+    The returned value for each time is either the monthly value or an
+    interpolated value to smooth the transition between months.
+    Interpolation is done on the day of year as determined by UTC.
     """
 
     # The .h5 file 'LinkeTurbidities.h5' contains a single 2160 x 4320 x 12
@@ -200,7 +207,7 @@ def lookup_linke_turbidity(time, latitude, longitude, filepath=None,
     if interp_turbidity:
         linke_turbidity = _interpolate_turbidity(lts, time)
     else:
-        months = time.month - 1
+        months = tools._pandas_to_utc(time).month - 1
         linke_turbidity = pd.Series(lts[months], index=time)
 
     linke_turbidity /= 20.
@@ -246,14 +253,11 @@ def _interpolate_turbidity(lts, time):
     # Jan 1 - Jan 15 and Dec 16 - Dec 31.
     lts_concat = np.concatenate([[lts[-1]], lts, [lts[0]]])
 
-    # handle leap years
-    try:
-        isleap = time.is_leap_year
-    except AttributeError:
-        year = time.year
-        isleap = _is_leap_year(year)
+    time_utc = tools._pandas_to_utc(time)
+
+    isleap = time_utc.is_leap_year
 
-    dayofyear = time.dayofyear
+    dayofyear = time_utc.dayofyear
     days_leap = _calendar_month_middles(2016)
     days_no_leap = _calendar_month_middles(2015)
 

pvlib/solarposition.py

@@ -23,7 +23,7 @@
 import scipy.optimize as so
 import warnings
 
-from pvlib import atmosphere
+from pvlib import atmosphere, tools
 from pvlib.tools import datetime_to_djd, djd_to_datetime
 
 
@@ -199,11 +199,7 @@ def spa_c(time, latitude, longitude, pressure=101325, altitude=0,
         raise ImportError('Could not import built-in SPA calculator. ' +
                           'You may need to recompile the SPA code.')
 
-    # if localized, convert to UTC. otherwise, assume UTC.
-    try:
-        time_utc = time.tz_convert('UTC')
-    except TypeError:
-        time_utc = time
+    time_utc = tools._pandas_to_utc(time)
 
     spa_out = []
 
@@ -378,7 +374,9 @@ def spa_python(time, latitude, longitude,
 
     spa = _spa_python_import(how)
 
-    delta_t = delta_t or spa.calculate_deltat(time.year, time.month)
+    if not delta_t:
+        time_utc = tools._pandas_to_utc(time)
+        delta_t = spa.calculate_deltat(time_utc.year, time_utc.month)
 
     app_zenith, zenith, app_elevation, elevation, azimuth, eot = \
         spa.solar_position(unixtime, lat, lon, elev, pressure, temperature,
@@ -452,12 +450,13 @@ def sun_rise_set_transit_spa(times, latitude, longitude, how='numpy',
         raise ValueError('times must be localized')
 
     # must convert to midnight UTC on day of interest
-    utcday = pd.DatetimeIndex(times.date).tz_localize('UTC')
-    unixtime = _datetime_to_unixtime(utcday)
+    times_utc = times.tz_convert('UTC')
+    unixtime = _datetime_to_unixtime(times_utc.normalize())
 
     spa = _spa_python_import(how)
 
-    delta_t = delta_t or spa.calculate_deltat(times.year, times.month)
+    if not delta_t:
+        delta_t = spa.calculate_deltat(times_utc.year, times_utc.month)
 
     transit, sunrise, sunset = spa.transit_sunrise_sunset(
         unixtime, lat, lon, delta_t, numthreads)
@@ -581,12 +580,11 @@ def sun_rise_set_transit_ephem(times, latitude, longitude,
     sunrise = []
     sunset = []
     trans = []
-    for thetime in times:
-        thetime = thetime.to_pydatetime()
+    for thetime in tools._pandas_to_utc(times):
         # older versions of pyephem ignore timezone when converting to its
         # internal datetime format, so convert to UTC here to support
         # all versions.  GH #1449
-        obs.date = ephem.Date(thetime.astimezone(dt.timezone.utc))
+        obs.date = ephem.Date(thetime)
         sunrise.append(_ephem_to_timezone(rising(sun), tzinfo))
         sunset.append(_ephem_to_timezone(setting(sun), tzinfo))
         trans.append(_ephem_to_timezone(transit(sun), tzinfo))
@@ -644,11 +642,7 @@ def pyephem(time, latitude, longitude, altitude=0, pressure=101325,
     except ImportError:
         raise ImportError('PyEphem must be installed')
 
-    # if localized, convert to UTC. otherwise, assume UTC.
-    try:
-        time_utc = time.tz_convert('UTC')
-    except TypeError:
-        time_utc = time
+    time_utc = tools._pandas_to_utc(time)
 
     sun_coords = pd.DataFrame(index=time)
 
@@ -765,11 +759,7 @@ def ephemeris(time, latitude, longitude, pressure=101325, temperature=12):
     # the SPA algorithm needs time to be expressed in terms of
     # decimal UTC hours of the day of the year.
 
-    # if localized, convert to UTC. otherwise, assume UTC.
-    try:
-        time_utc = time.tz_convert('UTC')
-    except TypeError:
-        time_utc = time
+    time_utc = tools._pandas_to_utc(time)
 
     # strip out the day of the year and calculate the decimal hour
     DayOfYear = time_utc.dayofyear
@@ -956,7 +946,10 @@ def pyephem_earthsun_distance(time):
 
     sun = ephem.Sun()
     earthsun = []
-    for thetime in time:
+    for thetime in tools._pandas_to_utc(time):
+        # older versions of pyephem ignore timezone when converting to its
+        # internal datetime format, so convert to UTC here to support
+        # all versions.  GH #1449
         sun.compute(ephem.Date(thetime))
         earthsun.append(sun.earth_distance)
 
@@ -1013,7 +1006,9 @@ def nrel_earthsun_distance(time, how='numpy', delta_t=67.0, numthreads=4):
 
     spa = _spa_python_import(how)
 
-    delta_t = delta_t or spa.calculate_deltat(time.year, time.month)
+    if not delta_t:
+        time_utc = tools._pandas_to_utc(time)
+        delta_t = spa.calculate_deltat(time_utc.year, time_utc.month)
 
     dist = spa.earthsun_distance(unixtime, delta_t, numthreads)
 
@@ -1386,22 +1381,26 @@ def hour_angle(times, longitude, equation_of_time):
     equation_of_time_spencer71
     equation_of_time_pvcdrom
     """
+
+    # times must be localized
+    if not times.tz:
+        raise ValueError('times must be localized')
+
     # hours - timezone = (times - normalized_times) - (naive_times - times)
-    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
+    hrs_minus_tzs = _times_to_hours_after_local_midnight(times) - 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.)
+    return 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"""
-    if times.tz is None:
-        times = times.tz_localize('utc')
+
+    # times must be localized
+    if not times.tz:
+        raise ValueError('times must be localized')
+
     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)
@@ -1411,18 +1410,26 @@ 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
+
+    # times must be localized
+    if not times.tz:
+        raise ValueError('times must be localized')
+
+    # normalize local times to previous local midnight and add the hours until
     # sunrise, sunset, and transit
-    return pd.DatetimeIndex(
-        naive_times.normalize() + pd.to_timedelta(hours, unit='h'), tz=tz_info)
+    return times.normalize() + pd.to_timedelta(hours, unit='h')
 
 
 def _times_to_hours_after_local_midnight(times):
     """convert local pandas datetime indices to array of hours as floats"""
-    times = times.tz_localize(None)
+
+    # times must be localized
+    if not times.tz:
+        raise ValueError('times must be localized')
+
     hrs = (times - times.normalize()) / pd.Timedelta('1h')
+
+    # ensure array return instead of a version-dependent pandas <T>Index
     return np.array(hrs)
 
 
@@ -1468,6 +1475,11 @@ def sun_rise_set_transit_geometric(times, latitude, longitude, declination,
        CRC Press (2012)
 
     """
+
+    # times must be localized
+    if not times.tz:
+        raise ValueError('times must be localized')
+
     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

pvlib/spectrum/spectrl2.py

@@ -288,7 +288,7 @@ def spectrl2(apparent_zenith, aoi, surface_tilt, ground_albedo,
                 aerosol_turbidity_500nm, scattering_albedo_400nm, alpha,
                 wavelength_variation_factor, aerosol_asymmetry_factor]))
 
-        dayofyear = original_index.dayofyear.values
+        dayofyear = pvlib.tools._pandas_to_doy(original_index).values
 
     if not is_pandas and dayofyear is None:
         raise ValueError('dayofyear must be specified if not using pandas '