diff --git a/docs/examples/irradiance-transposition/plot_seasonal_tilt.py b/docs/examples/irradiance-transposition/plot_seasonal_tilt.py index afe610f6d2..dc4b433412 100644 --- a/docs/examples/irradiance-transposition/plot_seasonal_tilt.py +++ b/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, diff --git a/docs/sphinx/source/whatsnew/v0.11.1.rst b/docs/sphinx/source/whatsnew/v0.11.1.rst index 7c8f01743e..7734cf3006 100644 --- a/docs/sphinx/source/whatsnew/v0.11.1.rst +++ b/docs/sphinx/source/whatsnew/v0.11.1.rst @@ -33,7 +33,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 ~~~~~~~ @@ -61,6 +73,7 @@ Requirements Contributors ~~~~~~~~~~~~ * Echedey Luis (:ghuser:`echedey-ls`) +* Yunho Kee (:ghuser:`yhkee0404`) * Chris Deline (:ghuser:`cdeline`) * Ioannis Sifnaios (:ghuser:`IoannisSifnaios`) * Leonardo Micheli (:ghuser:`lmicheli`) diff --git a/docs/tutorials/irradiance.ipynb b/docs/tutorials/irradiance.ipynb index 09b1323658..251d8b5eb8 100644 --- a/docs/tutorials/irradiance.ipynb +++ b/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", diff --git a/pvlib/clearsky.py b/pvlib/clearsky.py index 62d83bc1f4..626e4080ef 100644 --- a/pvlib/clearsky.py +++ b/pvlib/clearsky.py @@ -169,6 +169,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 @@ -201,7 +208,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. @@ -247,14 +254,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) diff --git a/pvlib/solarposition.py b/pvlib/solarposition.py index 4b76204063..5a3cec2dfc 100644 --- a/pvlib/solarposition.py +++ b/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 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 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 diff --git a/pvlib/spectrum/spectrl2.py b/pvlib/spectrum/spectrl2.py index ef25cb3787..b022d79371 100644 --- a/pvlib/spectrum/spectrl2.py +++ b/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 ' diff --git a/pvlib/tests/test_clearsky.py b/pvlib/tests/test_clearsky.py index 4e353e997c..509db83469 100644 --- a/pvlib/tests/test_clearsky.py +++ b/pvlib/tests/test_clearsky.py @@ -749,6 +749,7 @@ def test_bird(): tz = -7 # test timezone gmt_tz = pytz.timezone('Etc/GMT%+d' % -(tz)) times = times.tz_localize(gmt_tz) # set timezone + times_utc = times.tz_convert('UTC') # match test data from BIRD_08_16_2012.xls latitude = 40. longitude = -105. @@ -759,9 +760,9 @@ def test_bird(): aod_380nm = 0.15 b_a = 0.85 alb = 0.2 - eot = solarposition.equation_of_time_spencer71(times.dayofyear) + eot = solarposition.equation_of_time_spencer71(times_utc.dayofyear) hour_angle = solarposition.hour_angle(times, longitude, eot) - 0.5 * 15. - declination = solarposition.declination_spencer71(times.dayofyear) + declination = solarposition.declination_spencer71(times_utc.dayofyear) zenith = solarposition.solar_zenith_analytical( np.deg2rad(latitude), np.deg2rad(hour_angle), declination ) @@ -777,32 +778,35 @@ def test_bird(): Eb, Ebh, Gh, Dh = (irrads[_] for _ in field_names) data_path = DATA_DIR / 'BIRD_08_16_2012.csv' testdata = pd.read_csv(data_path, usecols=range(1, 26), header=1).dropna() - testdata.index = times[1:48] - assert np.allclose(testdata['DEC'], np.rad2deg(declination[1:48])) - assert np.allclose(testdata['EQT'], eot[1:48], rtol=1e-4) - assert np.allclose(testdata['Hour Angle'], hour_angle[1:48]) - assert np.allclose(testdata['Zenith Ang'], zenith[1:48]) + testdata[['DEC', 'EQT']] = testdata[['DEC', 'EQT']].shift(tz) + testdata = testdata[:tz] + end = 48 + tz + testdata.index = times[1:end] + assert np.allclose(testdata['DEC'], np.rad2deg(declination[1:end])) + assert np.allclose(testdata['EQT'], eot[1:end], rtol=1e-4) + assert np.allclose(testdata['Hour Angle'], hour_angle[1:end], rtol=1e-2) + assert np.allclose(testdata['Zenith Ang'], zenith[1:end], rtol=1e-2) dawn = zenith < 88. dusk = testdata['Zenith Ang'] < 88. am = pd.Series(np.where(dawn, airmass, 0.), index=times).fillna(0.0) assert np.allclose( - testdata['Air Mass'].where(dusk, 0.), am[1:48], rtol=1e-3 + testdata['Air Mass'].where(dusk, 0.), am[1:end], rtol=1e-3 ) direct_beam = pd.Series(np.where(dawn, Eb, 0.), index=times).fillna(0.) assert np.allclose( - testdata['Direct Beam'].where(dusk, 0.), direct_beam[1:48], rtol=1e-3 + testdata['Direct Beam'].where(dusk, 0.), direct_beam[1:end], rtol=1e-3 ) direct_horz = pd.Series(np.where(dawn, Ebh, 0.), index=times).fillna(0.) assert np.allclose( - testdata['Direct Hz'].where(dusk, 0.), direct_horz[1:48], rtol=1e-3 + testdata['Direct Hz'].where(dusk, 0.), direct_horz[1:end], rtol=1e-3 ) global_horz = pd.Series(np.where(dawn, Gh, 0.), index=times).fillna(0.) assert np.allclose( - testdata['Global Hz'].where(dusk, 0.), global_horz[1:48], rtol=1e-3 + testdata['Global Hz'].where(dusk, 0.), global_horz[1:end], rtol=1e-3 ) diffuse_horz = pd.Series(np.where(dawn, Dh, 0.), index=times).fillna(0.) assert np.allclose( - testdata['Dif Hz'].where(dusk, 0.), diffuse_horz[1:48], rtol=1e-3 + testdata['Dif Hz'].where(dusk, 0.), diffuse_horz[1:end], rtol=1e-3 ) # repeat test with albedo as a Series alb_series = pd.Series(0.2, index=times) @@ -813,19 +817,19 @@ def test_bird(): Eb, Ebh, Gh, Dh = (irrads[_] for _ in field_names) direct_beam = pd.Series(np.where(dawn, Eb, 0.), index=times).fillna(0.) assert np.allclose( - testdata['Direct Beam'].where(dusk, 0.), direct_beam[1:48], rtol=1e-3 + testdata['Direct Beam'].where(dusk, 0.), direct_beam[1:end], rtol=1e-3 ) direct_horz = pd.Series(np.where(dawn, Ebh, 0.), index=times).fillna(0.) assert np.allclose( - testdata['Direct Hz'].where(dusk, 0.), direct_horz[1:48], rtol=1e-3 + testdata['Direct Hz'].where(dusk, 0.), direct_horz[1:end], rtol=1e-3 ) global_horz = pd.Series(np.where(dawn, Gh, 0.), index=times).fillna(0.) assert np.allclose( - testdata['Global Hz'].where(dusk, 0.), global_horz[1:48], rtol=1e-3 + testdata['Global Hz'].where(dusk, 0.), global_horz[1:end], rtol=1e-3 ) diffuse_horz = pd.Series(np.where(dawn, Dh, 0.), index=times).fillna(0.) assert np.allclose( - testdata['Dif Hz'].where(dusk, 0.), diffuse_horz[1:48], rtol=1e-3 + testdata['Dif Hz'].where(dusk, 0.), diffuse_horz[1:end], rtol=1e-3 ) # test keyword parameters @@ -835,22 +839,25 @@ def test_bird(): Eb2, Ebh2, Gh2, Dh2 = (irrads2[_] for _ in field_names) data_path = DATA_DIR / 'BIRD_08_16_2012_patm.csv' testdata2 = pd.read_csv(data_path, usecols=range(1, 26), header=1).dropna() - testdata2.index = times[1:48] + testdata2[['DEC', 'EQT']] = testdata2[['DEC', 'EQT']].shift(tz) + testdata2 = testdata2[:tz] + testdata2.index = times[1:end] direct_beam2 = pd.Series(np.where(dawn, Eb2, 0.), index=times).fillna(0.) assert np.allclose( - testdata2['Direct Beam'].where(dusk, 0.), direct_beam2[1:48], rtol=1e-3 + testdata2['Direct Beam'].where(dusk, 0.), direct_beam2[1:end], + rtol=1e-3 ) direct_horz2 = pd.Series(np.where(dawn, Ebh2, 0.), index=times).fillna(0.) assert np.allclose( - testdata2['Direct Hz'].where(dusk, 0.), direct_horz2[1:48], rtol=1e-3 + testdata2['Direct Hz'].where(dusk, 0.), direct_horz2[1:end], rtol=1e-3 ) global_horz2 = pd.Series(np.where(dawn, Gh2, 0.), index=times).fillna(0.) assert np.allclose( - testdata2['Global Hz'].where(dusk, 0.), global_horz2[1:48], rtol=1e-3 + testdata2['Global Hz'].where(dusk, 0.), global_horz2[1:end], rtol=1e-3 ) diffuse_horz2 = pd.Series(np.where(dawn, Dh2, 0.), index=times).fillna(0.) assert np.allclose( - testdata2['Dif Hz'].where(dusk, 0.), diffuse_horz2[1:48], rtol=1e-3 + testdata2['Dif Hz'].where(dusk, 0.), diffuse_horz2[1:end], rtol=1e-3 ) # test scalars just at noon # XXX: calculations start at 12am so noon is at index = 12 diff --git a/pvlib/tests/test_irradiance.py b/pvlib/tests/test_irradiance.py index 0eb951f9a1..3636b3b400 100644 --- a/pvlib/tests/test_irradiance.py +++ b/pvlib/tests/test_irradiance.py @@ -603,11 +603,11 @@ def test_poa_components(irrad_data, ephem_data, dni_et, relative_airmass): @pytest.mark.parametrize('pressure,expected', [ (93193, [[830.46567, 0.79742, 0.93505], - [676.09497, 0.63776, 3.02102]]), + [676.18340, 0.63782, 3.02102]]), (None, [[868.72425, 0.79742, 1.01664], - [680.66679, 0.63776, 3.28463]]), + [680.73800, 0.63782, 3.28463]]), (101325, [[868.72425, 0.79742, 1.01664], - [680.66679, 0.63776, 3.28463]]) + [680.73800, 0.63782, 3.28463]]) ]) def test_disc_value(pressure, expected): # see GH 449 for pressure=None vs. 101325. @@ -1080,7 +1080,7 @@ def test_dirindex(times): pressure=pressure, use_delta_kt_prime=True, temp_dew=tdew).values - expected_out = np.array([np.nan, 0., 748.31562753, 630.72592644]) + expected_out = np.array([np.nan, 0., 748.31562800, 630.73752100]) tolerance = 1e-8 assert np.allclose(out, expected_out, rtol=tolerance, atol=0, diff --git a/pvlib/tests/test_solarposition.py b/pvlib/tests/test_solarposition.py index 472383acce..6fae71ef79 100644 --- a/pvlib/tests/test_solarposition.py +++ b/pvlib/tests/test_solarposition.py @@ -139,7 +139,8 @@ def test_spa_python_numpy_physical_dst(expected_solpos, golden): assert_frame_equal(expected_solpos, ephem_data[expected_solpos.columns]) -def test_sun_rise_set_transit_spa(expected_rise_set_spa, golden): +@pytest.mark.parametrize('delta_t', [65.0, None]) +def test_sun_rise_set_transit_spa(expected_rise_set_spa, golden, delta_t): # solution from NREL SAP web calculator south = Location(-35.0, 0.0, tz='UTC') times = pd.DatetimeIndex([datetime.datetime(1996, 7, 5, 0), @@ -160,7 +161,7 @@ def test_sun_rise_set_transit_spa(expected_rise_set_spa, golden): result = solarposition.sun_rise_set_transit_spa(times, south.latitude, south.longitude, - delta_t=65.0) + delta_t=delta_t) result_rounded = pd.DataFrame(index=result.index) # need to iterate because to_datetime does not accept 2D data # the rounding fails on pandas < 0.17 @@ -172,7 +173,7 @@ def test_sun_rise_set_transit_spa(expected_rise_set_spa, golden): # test for Golden, CO compare to NREL SPA result = solarposition.sun_rise_set_transit_spa( expected_rise_set_spa.index, golden.latitude, golden.longitude, - delta_t=65.0) + delta_t=delta_t) # round to nearest minute result_rounded = pd.DataFrame(index=result.index) @@ -529,17 +530,18 @@ def test_get_solarposition_method_pyephem(expected_solpos, golden): assert_frame_equal(expected_solpos, ephem_data[expected_solpos.columns]) -def test_nrel_earthsun_distance(): +@pytest.mark.parametrize('delta_t', [64.0, None]) +def test_nrel_earthsun_distance(delta_t): times = pd.DatetimeIndex([datetime.datetime(2015, 1, 2), datetime.datetime(2015, 8, 2)] ).tz_localize('MST') - result = solarposition.nrel_earthsun_distance(times, delta_t=64.0) + result = solarposition.nrel_earthsun_distance(times, delta_t=delta_t) expected = pd.Series(np.array([0.983289204601, 1.01486146446]), index=times) assert_series_equal(expected, result) times = datetime.datetime(2015, 1, 2) - result = solarposition.nrel_earthsun_distance(times, delta_t=64.0) + result = solarposition.nrel_earthsun_distance(times, delta_t=delta_t) expected = pd.Series(np.array([0.983289204601]), index=pd.DatetimeIndex([times, ])) assert_series_equal(expected, result) @@ -579,19 +581,20 @@ def test_declination(): def test_analytical_zenith(): times = pd.date_range(start="1/1/2015 0:00", end="12/31/2015 23:00", freq="h").tz_localize('Etc/GMT+8') + times_utc = times.tz_convert('UTC') lat, lon = 37.8, -122.25 lat_rad = np.deg2rad(lat) output = solarposition.spa_python(times, lat, lon, 100) solar_zenith = np.deg2rad(output['zenith']) # spa # spencer - eot = solarposition.equation_of_time_spencer71(times.dayofyear) + eot = solarposition.equation_of_time_spencer71(times_utc.dayofyear) hour_angle = np.deg2rad(solarposition.hour_angle(times, lon, eot)) - decl = solarposition.declination_spencer71(times.dayofyear) + decl = solarposition.declination_spencer71(times_utc.dayofyear) zenith_1 = solarposition.solar_zenith_analytical(lat_rad, hour_angle, decl) # pvcdrom and cooper - eot = solarposition.equation_of_time_pvcdrom(times.dayofyear) + eot = solarposition.equation_of_time_pvcdrom(times_utc.dayofyear) hour_angle = np.deg2rad(solarposition.hour_angle(times, lon, eot)) - decl = solarposition.declination_cooper69(times.dayofyear) + decl = solarposition.declination_cooper69(times_utc.dayofyear) zenith_2 = solarposition.solar_zenith_analytical(lat_rad, hour_angle, decl) assert np.allclose(zenith_1, solar_zenith, atol=0.015) assert np.allclose(zenith_2, solar_zenith, atol=0.025) @@ -600,22 +603,23 @@ def test_analytical_zenith(): def test_analytical_azimuth(): times = pd.date_range(start="1/1/2015 0:00", end="12/31/2015 23:00", freq="h").tz_localize('Etc/GMT+8') + times_utc = times.tz_convert('UTC') lat, lon = 37.8, -122.25 lat_rad = np.deg2rad(lat) output = solarposition.spa_python(times, lat, lon, 100) solar_azimuth = np.deg2rad(output['azimuth']) # spa solar_zenith = np.deg2rad(output['zenith']) # spencer - eot = solarposition.equation_of_time_spencer71(times.dayofyear) + eot = solarposition.equation_of_time_spencer71(times_utc.dayofyear) hour_angle = np.deg2rad(solarposition.hour_angle(times, lon, eot)) - decl = solarposition.declination_spencer71(times.dayofyear) + decl = solarposition.declination_spencer71(times_utc.dayofyear) zenith = solarposition.solar_zenith_analytical(lat_rad, hour_angle, decl) azimuth_1 = solarposition.solar_azimuth_analytical(lat_rad, hour_angle, decl, zenith) # pvcdrom and cooper - eot = solarposition.equation_of_time_pvcdrom(times.dayofyear) + eot = solarposition.equation_of_time_pvcdrom(times_utc.dayofyear) hour_angle = np.deg2rad(solarposition.hour_angle(times, lon, eot)) - decl = solarposition.declination_cooper69(times.dayofyear) + decl = solarposition.declination_cooper69(times_utc.dayofyear) zenith = solarposition.solar_zenith_analytical(lat_rad, hour_angle, decl) azimuth_2 = solarposition.solar_azimuth_analytical(lat_rad, hour_angle, decl, zenith) @@ -665,21 +669,45 @@ def test_hour_angle(): '2015-01-02 12:04:44.6340' ]).tz_localize('Etc/GMT+7') eot = np.array([-3.935172, -4.117227, -4.026295]) - hours = solarposition.hour_angle(times, longitude, eot) + hourangle = solarposition.hour_angle(times, longitude, eot) expected = (-70.682338, 70.72118825000001, 0.000801250) # FIXME: there are differences from expected NREL SPA calculator values # sunrise: 4 seconds, sunset: 48 seconds, transit: 0.2 seconds # but the differences may be due to other SPA input parameters - assert np.allclose(hours, expected) + assert np.allclose(hourangle, expected) + + hours = solarposition._hour_angle_to_hours( + times, hourangle, longitude, eot) + result = solarposition._times_to_hours_after_local_midnight(times) + assert np.allclose(result, hours) + + result = solarposition._local_times_from_hours_since_midnight(times, hours) + assert result.equals(times) + + times = times.tz_convert(None) + with pytest.raises(ValueError): + solarposition.hour_angle(times, longitude, eot) + with pytest.raises(ValueError): + solarposition._hour_angle_to_hours(times, hourangle, longitude, eot) + with pytest.raises(ValueError): + solarposition._times_to_hours_after_local_midnight(times) + with pytest.raises(ValueError): + solarposition._local_times_from_hours_since_midnight(times, hours) def test_sun_rise_set_transit_geometric(expected_rise_set_spa, golden_mst): """Test geometric calculations for sunrise, sunset, and transit times""" times = expected_rise_set_spa.index + times_utc = times.tz_convert('UTC') latitude = golden_mst.latitude longitude = golden_mst.longitude - eot = solarposition.equation_of_time_spencer71(times.dayofyear) # minutes - decl = solarposition.declination_spencer71(times.dayofyear) # radians + eot = solarposition.equation_of_time_spencer71( + times_utc.dayofyear) # minutes + decl = solarposition.declination_spencer71(times_utc.dayofyear) # radians + with pytest.raises(ValueError): + solarposition.sun_rise_set_transit_geometric( + times.tz_convert(None), latitude=latitude, longitude=longitude, + declination=decl, equation_of_time=eot) sr, ss, st = solarposition.sun_rise_set_transit_geometric( times, latitude=latitude, longitude=longitude, declination=decl, equation_of_time=eot) @@ -781,7 +809,7 @@ def test_nrel_earthsun_distance_microsecond_index(tz): @requires_pandas_2_0 -@pytest.mark.parametrize('tz', [None, 'utc', 'US/Eastern']) +@pytest.mark.parametrize('tz', ['utc', 'US/Eastern']) def test_hour_angle_microsecond_index(tz): # https://github.com/pvlib/pvlib-python/issues/1932 @@ -813,7 +841,7 @@ def test_rise_set_transit_spa_microsecond_index(tz): @requires_pandas_2_0 -@pytest.mark.parametrize('tz', [None, 'utc', 'US/Eastern']) +@pytest.mark.parametrize('tz', ['utc', 'US/Eastern']) def test_rise_set_transit_geometric_microsecond_index(tz): # https://github.com/pvlib/pvlib-python/issues/1932 diff --git a/pvlib/tools.py b/pvlib/tools.py index 3d766b6f72..c8d4d6e309 100644 --- a/pvlib/tools.py +++ b/pvlib/tools.py @@ -206,8 +206,32 @@ def _pandas_to_doy(pd_object): Returns ------- dayofyear + + Notes + ----- + Day of year is determined using UTC, since pandas uses local hour + """ + return _pandas_to_utc(pd_object).dayofyear + + +def _pandas_to_utc(pd_object): """ - return pd_object.dayofyear + Converts a pandas datetime-like object to UTC, if localized. + Otherwise, assume UTC. + + Parameters + ---------- + pd_object : DatetimeIndex or Timestamp + + Returns + ------- + pandas object localized to or assumed to be UTC. + """ + try: + pd_object_utc = pd_object.tz_convert('UTC') + except TypeError: + pd_object_utc = pd_object + return pd_object_utc def _doy_to_datetimeindex(doy, epoch_year=2014): @@ -230,7 +254,7 @@ def _doy_to_datetimeindex(doy, epoch_year=2014): def _datetimelike_scalar_to_doy(time): - return pd.DatetimeIndex([pd.Timestamp(time)]).dayofyear + return _pandas_to_doy(_datetimelike_scalar_to_datetimeindex(time)) def _datetimelike_scalar_to_datetimeindex(time):