diff --git a/docs/sphinx/source/reference/solarposition.rst b/docs/sphinx/source/reference/solarposition.rst index dbe93d5392..555151ef7b 100644 --- a/docs/sphinx/source/reference/solarposition.rst +++ b/docs/sphinx/source/reference/solarposition.rst @@ -17,6 +17,7 @@ unless you know that you need a different function. solarposition.get_solarposition solarposition.spa_python solarposition.ephemeris + solarposition.geometric solarposition.pyephem solarposition.spa_c diff --git a/docs/sphinx/source/whatsnew/v0.16.2.rst b/docs/sphinx/source/whatsnew/v0.16.2.rst index 64132a2547..3f335282b6 100644 --- a/docs/sphinx/source/whatsnew/v0.16.2.rst +++ b/docs/sphinx/source/whatsnew/v0.16.2.rst @@ -27,6 +27,10 @@ Enhancements * Add ``dataset`` parameter to :py:func:`~pvlib.iotools.get_era5` allowing for choosing between ERA5 and ERA5-Land datasets. (:pull:`2780`) +* Add :py:func:`pvlib.solarposition.geometric`, a fast but approximate solar + position calculation based on spherical geometry, and make it available as + ``method='geometric'`` in :py:func:`pvlib.solarposition.get_solarposition`. + (:issue:`2914`, :pull:`2915`) Documentation ~~~~~~~~~~~~~ @@ -69,3 +73,4 @@ Contributors * jason-writer (:ghuser:`jason-writer`) * :ghuser:`bfs2021` * :ghuser:`xiajiahang` +* Divyanshi Singh (:ghuser:`divyanshisingh987456321`) diff --git a/pvlib/solarposition.py b/pvlib/solarposition.py index 4441ed01f1..61c9568f8a 100644 --- a/pvlib/solarposition.py +++ b/pvlib/solarposition.py @@ -69,6 +69,9 @@ def get_solarposition(time, latitude, longitude, 'nrel_c' uses the NLR SPA C code [3]: :py:func:`spa_c` + 'geometric' uses a fast but approximate spherical-geometry + calculation: :py:func:`geometric` + temperature : float, default 12 Degrees C. @@ -119,6 +122,8 @@ def get_solarposition(time, latitude, longitude, elif method == 'ephemeris': ephem_df = ephemeris(time, latitude, longitude, pressure, temperature, **kwargs) + elif method == 'geometric': + ephem_df = geometric(time, latitude, longitude, **kwargs) else: raise ValueError('Invalid solar position method') @@ -878,6 +883,83 @@ def ephemeris(time, latitude, longitude, pressure=101325.0, temperature=12.0): return result +def geometric(time, latitude, longitude): + """ + Fast, approximate solar position from spherical geometry. + + The solar declination and the equation of time are evaluated with the + Spencer (1971) Fourier series, and the zenith and azimuth angles follow + from spherical trigonometry. This is much faster than, but less + accurate than, :py:func:`spa_python`; use it only when speed matters more + than accuracy. + + .. warning:: + Atmospheric refraction is neglected, so ``apparent_zenith`` and + ``apparent_elevation`` are equal to ``zenith`` and ``elevation``. + + Parameters + ---------- + time : pandas.DatetimeIndex + Must be localized or UTC will be assumed. + latitude : float + Latitude in decimal degrees. Positive north of equator, negative + to south. + longitude : float + Longitude in decimal degrees. Positive east of prime meridian, + negative to west. + + Returns + ------- + DataFrame with the following columns: + + * apparent_elevation : equal to ``elevation`` (no refraction). + * elevation : elevation of the sun in decimal degrees, 0 = on + horizon. The complement of the zenith angle. + * azimuth : Azimuth of the sun in decimal degrees East of North. + * apparent_zenith : equal to ``zenith`` (no refraction). + * zenith : Solar zenith angle in decimal degrees. + + See Also + -------- + spa_python, ephemeris, declination_spencer71, + equation_of_time_spencer71, solar_zenith_analytical, + solar_azimuth_analytical + """ + time_utc = tools._pandas_to_utc(time) # naive times assumed to be UTC + + dayofyear = time_utc.dayofyear + declination = declination_spencer71(dayofyear) + eot = equation_of_time_spencer71(dayofyear) + utc_hours = (time_utc.hour + time_utc.minute / 60. + + time_utc.second / 3600. + time_utc.microsecond / 3.6e9) + # hour angle in degrees, as in hour_angle() but computed from UTC time + # (no timezone offset to apply) + hourangle = np.radians(15. * (np.asarray(utc_hours) - 12.) + longitude + + eot / 4.) + # wrap to [-pi, pi): the sign of the hour angle selects the azimuth branch + hourangle = np.mod(hourangle + np.pi, 2 * np.pi) - np.pi + latitude_rad = np.radians(latitude) + + zenith = solar_zenith_analytical(latitude_rad, hourangle, declination) + azimuth = solar_azimuth_analytical( + latitude_rad, hourangle, declination, zenith) + + zenith = np.degrees(zenith) + azimuth = np.mod(np.degrees(azimuth), 360) + elevation = 90 - zenith + + return pd.DataFrame( + { + "apparent_elevation": elevation, + "elevation": elevation, + "azimuth": azimuth, + "apparent_zenith": zenith, + "zenith": zenith, + }, + index=time, + ) + + def calc_time(lower_bound, upper_bound, latitude, longitude, attribute, value, altitude=0.0, pressure=101325.0, temperature=12.0, horizon='+0:00', xtol=1.0e-12): diff --git a/tests/test_solarposition.py b/tests/test_solarposition.py index 9dc0488fe5..fc04050c82 100644 --- a/tests/test_solarposition.py +++ b/tests/test_solarposition.py @@ -463,6 +463,58 @@ def test_ephemeris_physical_no_tz(expected_solpos, golden_mst): assert_frame_equal(expected_solpos, ephem_data[expected_solpos.columns]) +@pytest.mark.parametrize('latitude, longitude, tz', [ + (39.742476, -105.1786, 'America/Denver'), # west of Greenwich + (28.46, 77.03, 'Asia/Kolkata'), # east of Greenwich + (-33.87, 151.21, 'Australia/Sydney'), # southern hemisphere + (69.65, 18.96, 'Europe/Oslo'), # high latitude +]) +def test_geometric_vs_spa(latitude, longitude, tz): + # the geometric method is approximate: compare with SPA across a whole + # year, including evenings at western longitudes where the hour angle + # leaves [-180, 180] degrees before wrapping + times = pd.date_range('2026-01-01', '2026-12-31 23:00', freq='1h', tz=tz) + expected = solarposition.spa_python(times, latitude, longitude) + result = solarposition.geometric(times, latitude, longitude) + assert_allclose(result['zenith'], expected['zenith'], atol=1) + assert_allclose(result['elevation'], expected['elevation'], atol=1) + daytime = expected['elevation'] > 0 + azimuth_error = (result['azimuth'] - expected['azimuth'] + 180) % 360 - 180 + assert (azimuth_error[daytime].abs() < 2).all() + + +def test_geometric_output(golden): + result = solarposition.geometric(times_localized, golden.latitude, + golden.longitude) + assert result.index.equals(times_localized) + assert list(result.columns) == ['apparent_elevation', 'elevation', + 'azimuth', 'apparent_zenith', 'zenith'] + assert not result.isna().any().any() + assert result['azimuth'].between(0, 360).all() + assert_allclose(result['zenith'] + result['elevation'], 90) + # refraction is neglected + assert_series_equal(result['apparent_zenith'], result['zenith'], + check_names=False) + + +def test_geometric_no_tz(): + # times without a timezone are assumed to be UTC + times_naive = pd.date_range('2026-06-21 00:00', periods=24, freq='1h') + expected = solarposition.geometric(times_naive.tz_localize('UTC'), + 32.2, -111) + result = solarposition.geometric(times_naive, 32.2, -111) + assert_allclose(result.to_numpy(), expected.to_numpy()) + + +def test_get_solarposition_method_geometric(golden): + expected = solarposition.geometric(times_localized, golden.latitude, + golden.longitude) + result = solarposition.get_solarposition( + times_localized, golden.latitude, golden.longitude, + method='geometric') + assert_frame_equal(result, expected) + + def test_get_solarposition_error(golden): times = pd.date_range(datetime.datetime(2003, 10, 17, 13, 30, 30), periods=1, freq='D', tz=golden.tz)