diff --git a/docs/source/documentation.rst b/docs/source/documentation.rst index 9da7647..ec2b25d 100644 --- a/docs/source/documentation.rst +++ b/docs/source/documentation.rst @@ -11,6 +11,7 @@ Code documentation plotting.irradiance_colormap_and_norm plotting.plot_google_maps plotting.plot_intraday_heatmap + instrument.shadowband_correction_factor plotting.plot_shading_heatmap horizon.get_horizon_mines quality.bsrn_limits diff --git a/pyproject.toml b/pyproject.toml index e28dcbe..a3c5e27 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -47,8 +47,8 @@ dynamic = ["version"] [project.optional-dependencies] test = ["pytest>=7", "pytest-cov", "packaging"] doc = [ - "sphinx==9.1.0", - "myst-nb==1.4.0", + "sphinx==9.0.4", + "myst-nb==1.3.0", "sphinx-book-theme==1.1.4", "sphinx-copybutton==0.5.2", "pvlib==0.14.0", diff --git a/src/solarpy/__init__.py b/src/solarpy/__init__.py index d0c85fe..38c3819 100644 --- a/src/solarpy/__init__.py +++ b/src/solarpy/__init__.py @@ -11,5 +11,6 @@ quality, horizon, example, + instrument, iotools, ) diff --git a/src/solarpy/instrument/__init__.py b/src/solarpy/instrument/__init__.py new file mode 100644 index 0000000..5adce80 --- /dev/null +++ b/src/solarpy/instrument/__init__.py @@ -0,0 +1 @@ +from solarpy.instrument.shadowband import shadowband_correction_factor # noqa: F401 diff --git a/src/solarpy/instrument/shadowband.py b/src/solarpy/instrument/shadowband.py new file mode 100644 index 0000000..5eb2b38 --- /dev/null +++ b/src/solarpy/instrument/shadowband.py @@ -0,0 +1,128 @@ +""" +shadowband_correction.py +------------------------ +Shadowband diffuse-irradiance correction factor for the Kipp & Zonen CM 121 +shadow-ring pyranometer accessory. + +Reference +--------- +Kipp & Zonen *CM 121 B/C Shadow Ring Instruction Manual*, Appendix §6.1. +https://www.kippzonen.com/Download/46/CM121-B-C-Shadow-Ring-Manual +""" + +import numpy as np +import pvlib as pv +import pandas as pd +import datetime + + +def shadowband_correction_factor( + input_date, + latitude: float, + V: float = 0.185, +): + """Calculate the daily shadowband correction factor for a Kipp & Zonen CM 121. + + A shadow ring blocks a strip of the sky dome to prevent direct beam + radiation from reaching the pyranometer, so that only diffuse radiation + is measured. Because the ring also blocks some of the diffuse sky, a + correction factor C > 1 must be applied to recover the true diffuse + irradiance:: + + G_d_true = C * G_d_measured + + The factor depends on the ring geometry (angular half-width *V*), the + site latitude, and the solar declination on the measurement day. It is + derived analytically by integrating the fraction of the isotropic sky + hemisphere obscured by the ring over the daily arc of the sun. + + Parameters + ---------- + input_date : scalar or array-like + Date(s) for which to compute the factor. Accepts anything that + ``pandas.to_datetime`` understands: a ``pd.Timestamp``, a + ``datetime.date`` / ``datetime.datetime``, an ISO-8601 string, a + list/array of any of the above, or a ``pd.DatetimeIndex``. + Sub-daily time information is ignored; only the calendar date matters. + latitude : float + Observer latitude in decimal degrees. Positive = North, negative = South. + Valid range: −90 to +90. + V : float, optional + Angular width of the shadow band in radians. This is arc subtended by + the ring as seen from the sensor. Default is 0.185 rad (~10.6°), + which matches the standard Kipp & Zonen CM 121 ring geometry. + + Returns + ------- + float or array-like + Dimensionless correction factor C ≥ 1. + + Notes + ----- + The correction factor is computed as: + + .. math:: + + S = \\frac{2V \\cos\\delta}{\\pi} + \\bigl(U_0 \\sin\\phi\\sin\\delta + + \\sin U_0 \\cos\\phi\\cos\\delta\\bigr) + + C = \\frac{1}{1 - S} + + where + + * :math:`\\delta` – solar declination (Spencer 1971 approximation) + * :math:`\\phi` – site latitude + * :math:`U_0` – sunrise/sunset hour angle, + :math:`\\arccos(-\\tan\\phi\\,\\tan\\delta)`, clamped to [−1, 1] to + handle polar day/night conditions + """ + # ------------------------------------------------------------------ + # 1. Normalise input to pandas datetime so scalar and array paths + # share a single code route below. + # ------------------------------------------------------------------ + dates = pd.to_datetime(input_date) + + # ------------------------------------------------------------------ + # 2. Solar declination δ for each day of year. + # Spencer (1971) approximation, returned in radians by pvlib. + # ------------------------------------------------------------------ + day_of_year = dates.dayofyear + delta = pv.solarposition.declination_spencer71(day_of_year) # radians + + # ------------------------------------------------------------------ + # 3. Site latitude φ in radians. + # ------------------------------------------------------------------ + phi = np.deg2rad(latitude) + + # ------------------------------------------------------------------ + # 4. Sunrise/sunset hour angle U₀. + # cos U₀ = −tan φ · tan δ. + # Clamped to [−1, 1] so arccos is valid at the poles (midnight sun + # / polar night), where the tangent product can exceed ±1. + # ------------------------------------------------------------------ + cos_U0 = np.clip(-np.tan(phi) * np.tan(delta), -1.0, 1.0) + U0 = np.arccos(cos_U0) + + # ------------------------------------------------------------------ + # 5. Fraction S of the isotropic sky dome obscured by the ring. + # Derived by integrating the ring's shadow across the daily solar arc + # (Kipp & Zonen manual, Appendix §6.1). + # ------------------------------------------------------------------ + S = (2 * V * np.cos(delta) / np.pi) * ( + U0 * np.sin(phi) * np.sin(delta) + + np.sin(U0) * np.cos(phi) * np.cos(delta) + ) + + # ------------------------------------------------------------------ + # 6. Correction factor C = 1 / (1 − S). + # S is always < 1 for physically reasonable inputs, so C > 1. + # ------------------------------------------------------------------ + C = 1.0 / (1.0 - S) + + # ------------------------------------------------------------------ + # 7. Return a plain float for scalar input, Series for array input. + # ------------------------------------------------------------------ + if np.isscalar(input_date) or isinstance(input_date, (pd.Timestamp, datetime.date, datetime.datetime)): + return float(C) + return pd.Series(C, index=dates) diff --git a/tests/test_instrument.py b/tests/test_instrument.py new file mode 100644 index 0000000..7574e57 --- /dev/null +++ b/tests/test_instrument.py @@ -0,0 +1,58 @@ +import solarpy +import pandas as pd +import datetime +import pytest + + +def test_shadowband_correction(): + # Northern hemisphere, spring (near-zero declination) + assert solarpy.instrument.shadowband_correction_factor("2026-03-17", latitude=55.29007) == 1.0668908588370676 + # Northern hemisphere, spring - pd.Timestamp + assert solarpy.instrument.shadowband_correction_factor(pd.Timestamp("2026-03-17"), latitude=55.29007) == 1.0668908588370676 + # pd.Timestamp with time component - sub-daily info must be ignored + assert solarpy.instrument.shadowband_correction_factor(pd.Timestamp("2026-03-17 14:30:01"), latitude=55.29007) == pytest.approx(1.0668908588370676) + # datetime.date scalar + assert solarpy.instrument.shadowband_correction_factor(datetime.date(2026, 3, 17), latitude=55.29007) == pytest.approx(1.0668908588370676) + # datetime.datetime scalar + assert solarpy.instrument.shadowband_correction_factor(datetime.datetime(2026, 3, 17, 14, 30), latitude=55.29007) == pytest.approx(1.0668908588370676) + + # Northern hemisphere, summer + assert solarpy.instrument.shadowband_correction_factor("2026-06-21", latitude=55.29007) == pytest.approx(1.1408347415382885) + # Northern hemisphere, autumn + assert solarpy.instrument.shadowband_correction_factor("2026-10-21", latitude=55.29007) == pytest.approx(1.0418060380205798) + # Northern hemisphere, winter + assert solarpy.instrument.shadowband_correction_factor("2026-12-21", latitude=55.29007) == pytest.approx(1.0126112472695903) + + # Southern hemisphere + assert solarpy.instrument.shadowband_correction_factor("2026-10-21", latitude=-33.87) == pytest.approx(1.1282175055807726) + +def test_shadowband_correction_series(): + # Single value series + index = pd.to_datetime(["2026-03-17"]) + pd.testing.assert_series_equal( + solarpy.instrument.shadowband_correction_factor(index, latitude=55.29007), + pd.Series(1.0668908588370676, index=index)) + + # Multiple value series + index = pd.DatetimeIndex(["2026-03-17", "2026-10-21"]) + pd.testing.assert_series_equal( + solarpy.instrument.shadowband_correction_factor(index, latitude=55.29007), + pd.Series([1.0668908588370676, 1.0418060380205798], index=index)) + +def test_shadowband_correction_always_above_one(): + # Test if shadowband correction factor C > 1 + dates = pd.DatetimeIndex(["2026-03-20", "2026-06-21", "2026-09-22", "2026-12-21"]) + result = solarpy.instrument.shadowband_correction_factor(dates, latitude=55.29007) + assert (result > 1.0).all() + + +def test_shadowband_correction_v_zero(): + # With no ring there is nothing to correct for + assert solarpy.instrument.shadowband_correction_factor("2026-06-21", latitude=45.0, V=0.0) == pytest.approx(1.0) + +def test_shadowband_correction_polar(): + # Test that the clip works and it follows the table provided by Kipp & Zonen + # At the pole in winter the sun never rises — ring blocks nothing, C = 1 + assert solarpy.instrument.shadowband_correction_factor("2026-01-01", latitude=90.0) == pytest.approx(1.0) + # At the pole in summer the sun skims the horizon — C rises above 1 + assert solarpy.instrument.shadowband_correction_factor("2026-06-21", latitude=90.0) > 1.0