Fix cosine of solar zenith angle glitches (#1399)

### Pull Request Checklist: - [x] This PR addresses an already opened issue (for bug fixes / features) - This PR fixes #1110 - [x] Tests for the changes have been added (for bug fixes / features) - [x] (If applicable) Documentation has been added / updated (for bug fixes / features) - [x] CHANGES.rst has been updated (with summary of main changes) - [x] Link to issue (:issue:`number`) and pull request (:pull:`number`) has been added ### What kind of change does this PR introduce? Refactor of the function so we can use a ufunc, making it easier to implement the different cases. Earlier version was not handling the cyclical nature of angles correctly and had issues in the polar regions. ### Does this PR introduce a breaking change? Yes, the output has changed. I think it is now correct. The signature has also changed. ### Other information: PyWGBT was used as a comparison.
Ouranosinc · Jun 21, 2023 · 03d3376 · 03d3376
2 parents 3b2bff6 + c0a8244
commit 03d3376
Show file tree

Hide file tree

Showing 6 changed files with 257 additions and 156 deletions.
diff --git a/CHANGES.rst b/CHANGES.rst
@@ -35,6 +35,7 @@ Bug fixes
 * Fixed a bug in ``xclim.core.calendar.time_bnds`` when using ``DataArrayResample`` objects, caused by an upstream change in xarray 2023.5.0. (:issue:`1368`, :pull:`1377`).
 * `ensembles.change_significance` will returns NaNs when the input values are all NaNs, instead of failing. (:issue:`1379`, :pull:`1380`).
 * Accelerated import of xclim by caching the compilation of `guvectorize` functions. (:pull:`1378`).
+* Fixed many issues with ``xclim.indices.helpers.cosine_of_solar_zenith_angle``, the signature changed. (:issue:`1110`, :pull:`1399`).
 
 Internal changes
 ^^^^^^^^^^^^^^^^

diff --git a/tests/test_atmos.py b/tests/test_atmos.py
@@ -566,8 +566,8 @@ def test_mean_radiant_temperature(self, atmosds):
         rlus = dataset.rlus
 
         # Expected values
-        exp_sun = [np.nan, np.nan, np.nan, np.nan, np.nan]
-        exp_ins = [276.911, 274.742, 243.202, 268.012, 309.151]
+        exp_sun = [276.911, 274.742, 243.202, 268.012, 278.505]
+        exp_ins = [276.911, 274.742, 243.202, 268.012, 278.505]
         exp_avg = [276.911, 274.742, 243.202, 268.017, 278.512]
 
         mrt_sun = atmos.mean_radiant_temperature(rsds, rsus, rlds, rlus, stat="sunlit")

diff --git a/tests/test_helpers.py b/tests/test_helpers.py
@@ -5,7 +5,7 @@
 import pytest
 import xarray as xr
 
-from xclim.core.calendar import date_range, datetime_to_decimal_year
+from xclim.core.calendar import date_range
 from xclim.core.units import convert_units_to
 from xclim.indices import helpers
 
@@ -20,12 +20,8 @@ def test_solar_declinaton(method, rtol):
         dims=("time",),
     )
     exp = [-19.83, 20.64, 20.00]
-
-    day_angle = ((datetime_to_decimal_year(times) % 1) * 360).assign_attrs(
-        units="degree"
-    )
     np.testing.assert_allclose(
-        helpers.solar_declination(day_angle, method=method),
+        helpers.solar_declination(times, method=method),
         np.deg2rad(exp),
         atol=rtol * 2 * np.deg2rad(23.44),  # % of the possible range
     )
@@ -36,19 +32,17 @@ def test_extraterrestrial_radiation(method):
     # Expected values from https://www.engr.scu.edu/~emaurer/tools/calc_solar_cgi.pl
     # This source is not authoritative, thus the large rtol
     times = xr.DataArray(
-        pd.to_datetime(
-            ["1793-01-21T10:22:00", "1969-07-20T20:17:40", "2022-05-20T16:55:48"]
-        ),
+        xr.date_range("1900-01-01", "1900-01-03", freq="D"),
         dims=("time",),
         name="time",
     )
     lat = xr.DataArray(
-        [48.8656, 29.5519, 45.5435],
+        [48.8656, 29.5519, -54],
         dims=("time",),
         coords={"time": times},
         attrs={"units": "degree_north"},
     )
-    exp = [120.90, 477.51, 470.74]
+    exp = [99.06, 239.98, 520.01]
     np.testing.assert_allclose(
         convert_units_to(
             helpers.extraterrestrial_solar_radiation(times, lat, method=method), "W m-2"
@@ -91,3 +85,51 @@ def test_day_lengths(method):
                 np.testing.assert_allclose(
                     dl.sel(time=e[0]).transpose(), np.array(e[1]), rtol=2e-1
                 )
+
+
+def test_cosine_of_solar_zenith_angle():
+    time = xr.date_range("1900-01-01T00:30", "1900-01-03", freq="H")
+    time = xr.DataArray(time, dims=("time",), coords={"time": time}, name="time")
+    lat = xr.DataArray(
+        [0, 45, 70], dims=("site",), name="lat", attrs={"units": "degree_north"}
+    )
+    lon = xr.DataArray(
+        [-40, 0, 80], dims=("site",), name="lon", attrs={"units": "degree_east"}
+    )
+    dec = helpers.solar_declination(time)
+
+    czda = helpers.cosine_of_solar_zenith_angle(
+        time, dec, lat, lon, stat="average", sunlit=True
+    )
+    # Data Generated with PyWGBT
+    # raw = coszda(
+    #     (time + pd.Timedelta('30 m')).data,
+    #     convert_units_to(lat, 'rad').data[np.newaxis, :],
+    #     convert_units_to(lon, 'rad').data[np.newaxis, :],
+    #     1
+    # )
+    # exp_cza = xr.DataArray(raw, dims=('time', 'd', 'site'), coords={'lat': lat, 'lon': lon, 'time': time}).squeeze('d')
+    exp_czda = np.array(
+        [
+            [0.0, 0.0610457, 0.0],
+            [0.09999178, 0.18221077, 0.0],
+            [0.31387116, 0.285383, 0.0],
+            [0.52638271, 0.35026199, 0.0],
+            [0.70303168, 0.37242693, 0.0],
+        ]
+    )
+    np.testing.assert_allclose(czda[7:12, :], exp_czda, rtol=1e-3)
+
+    # Same code as above, but with function "cosza".
+    cza = helpers.cosine_of_solar_zenith_angle(
+        time, dec, lat, lon, stat="average", sunlit=False
+    )
+    exp_cza = np.array(
+        [
+            [-0.83153798, -0.90358335, -0.34065474],
+            [-0.90358299, -0.83874813, -0.26062708],
+            [-0.91405234, -0.73561867, -0.18790995],
+            [-0.86222963, -0.60121893, -0.12745608],
+        ]
+    )
+    np.testing.assert_allclose(cza[:4, :], exp_cza, rtol=1e-3)
diff --git a/tests/test_indices.py b/tests/test_indices.py
@@ -553,44 +553,44 @@ def test_standardized_precipitation_index(
                 1,
                 "gamma",
                 "APP",
-                [1.3750, 1.5776, 2.1033, -3.09, 0.8681],
+                [1.3750, 1.5776, 1.6806, -3.09, 0.8681],
                 2e-2,
             ),
             (
                 "MS",
                 12,
                 "gamma",
                 "APP",
-                [0.6242, 1.6479, 1.8005, 0.9857, 0.6706],
+                [0.6229, 1.6609, 1.8673, 1.0181, 0.6901],
                 2e-2,
             ),
-            ("MS", 1, "gamma", "ML", [1.38, 1.58, 2.1, -3.09, 0.868], 0.1),
+            ("MS", 1, "gamma", "ML", [1.38, 1.58, 2.1, -3.09, 0.868], 1e-1),
             ("MS", 1, "gamma", "ML", [1.4631, 1.6053, 2.1625, -3.09, 0.87996], 2e-2),
             (
                 "MS",
                 12,
                 "gamma",
                 "ML",
                 [0.62417, 1.6479, 1.8005, 0.98574, 0.67063],
-                0.06,
+                6e-2,
             ),
             (
                 "MS",
                 12,
                 "gamma",
                 "ML",
-                [0.65635, 1.5976, 1.7909, 0.98453, 0.66396],
+                [0.6511, 1.6146, 1.8580, 1.0140, 0.6878],
                 2e-2,
             ),
-            ("MS", 1, "fisk", "ML", [1.73, 1.51, 2.05, -3.09, 0.892], 0.35),
-            ("MS", 1, "fisk", "ML", [1.41675, 1.51174, 1.98248, -3.09, 0.942175], 2e-2),
-            ("MS", 12, "fisk", "ML", [0.71228, 1.5347, 1.6498, 1.0241, 0.72203], 0.03),
+            ("MS", 1, "fisk", "ML", [1.73, 1.51, 2.05, -3.09, 0.892], 3.5e-1),
+            ("MS", 1, "fisk", "ML", [1.4167, 1.5117, 2.0562, -3.09, 0.9422], 2e-2),
+            ("MS", 12, "fisk", "ML", [0.7041, 1.562985, 1.7041, 1.0388, 0.7165], 3e-2),
             (
                 "MS",
                 12,
                 "fisk",
                 "ML",
-                [0.70921, 1.55079, 1.6697, 1.02186, 0.695148],
+                [0.7041, 1.562985, 1.7041, 1.0388, 0.71645],
                 2e-2,
             ),
         ],
@@ -2982,7 +2982,9 @@ def test_baier_robertson(self, tasmin_series, tasmax_series, lat_series):
         tx = tasmax_series(np.array([10, 15, 20]) + 273.15).expand_dims(lat=lat)
 
         out = xci.potential_evapotranspiration(tn, tx, lat=lat, method="BR65")
-        np.testing.assert_allclose(out[0, 2], [3.861079 / 86400], rtol=1e-2)
+        np.testing.assert_allclose(
+            out.isel(lat=0, time=2), [3.861079 / 86400], rtol=1e-2
+        )
 
     def test_hargreaves(self, tasmin_series, tasmax_series, tas_series, lat_series):
         lat = lat_series([45])
@@ -2991,7 +2993,9 @@ def test_hargreaves(self, tasmin_series, tasmax_series, tas_series, lat_series):
         tm = tas_series(np.array([5, 10, 15]) + 273.15).expand_dims(lat=lat)
 
         out = xci.potential_evapotranspiration(tn, tx, tm, lat=lat, method="HG85")
-        np.testing.assert_allclose(out[0, 2], [3.962589 / 86400], rtol=1e-2)
+        np.testing.assert_allclose(
+            out.isel(lat=0, time=2), [3.962589 / 86400], rtol=1e-2
+        )
 
     def test_thornthwaite(self, tas_series, lat_series):
         lat = lat_series([45])
@@ -3007,15 +3011,19 @@ def test_thornthwaite(self, tas_series, lat_series):
 
         # find lat implicitly
         out = xci.potential_evapotranspiration(tas=tm, method="TW48")
-        np.testing.assert_allclose(out[0, 1], [42.7619242 / (86400 * 30)], rtol=1e-1)
+        np.testing.assert_allclose(
+            out.isel(lat=0, time=1), [42.7619242 / (86400 * 30)], rtol=1e-1
+        )
 
     def test_mcguinnessbordne(self, tasmin_series, tasmax_series, lat_series):
         lat = lat_series([45])
         tn = tasmin_series(np.array([0, 5, 10]) + 273.15).expand_dims(lat=lat)
         tx = tasmax_series(np.array([10, 15, 20]) + 273.15).expand_dims(lat=lat)
 
         out = xci.potential_evapotranspiration(tn, tx, lat=lat, method="MB05")
-        np.testing.assert_allclose(out[0, 2], [2.78253138816 / 86400], rtol=1e-2)
+        np.testing.assert_allclose(
+            out.isel(lat=0, time=2), [2.78253138816 / 86400], rtol=1e-2
+        )
 
     def test_allen(
         self,
@@ -3053,7 +3061,9 @@ def test_allen(
             sfcWind=sfcWind,
             method="FAO_PM98",
         )
-        np.testing.assert_allclose(out[0, 2], [1.208832768 / 86400], rtol=1e-2)
+        np.testing.assert_allclose(
+            out.isel(lat=0, time=2), [1.208832768 / 86400], rtol=1e-2
+        )
 
 
 def test_water_budget_from_tas(pr_series, tasmin_series, tasmax_series, lat_series):
@@ -3085,7 +3095,7 @@ def test_water_budget_from_tas(pr_series, tasmin_series, tasmax_series, lat_seri
 
     # find lat implicitly
     out = xci.water_budget(prm, tas=tm, method="TW48")
-    np.testing.assert_allclose(out[1, 0], [8.5746025 / 86400], rtol=2e-1)
+    np.testing.assert_allclose(out.isel(lat=0, time=1), [8.5746025 / 86400], rtol=2e-1)
 
 
 def test_water_budget(pr_series, evspsblpot_series):
@@ -3319,7 +3329,7 @@ def test_universal_thermal_climate_index(
 
 
 @pytest.mark.parametrize(
-    "stat,expected", [("sunlit", np.nan), ("instant", 295.3), ("average", 295.1)]
+    "stat,expected", [("sunlit", 295.0), ("instant", 294.9), ("average", 295.1)]
 )
 def test_mean_radiant_temperature(
     rsds_series,

diff --git a/xclim/indices/_conversion.py b/xclim/indices/_conversion.py
@@ -1359,8 +1359,8 @@ def potential_evapotranspiration(
     with :math:`a=0.0147` and :math:`b=0.07353`. The default parameters used here are calibrated for the UK,
     using the method described in :cite:t:`tanguy_historical_2018`.
 
-    Methods "BR65", "HG85" and "MB05" use an approximation of the extraterrestrial
-    radiation. See :py:func:`~xclim.indices._helpers.extraterrestrial_solar_radiation`.
+    Methods "BR65", "HG85" and "MB05" use an approximation of the extraterrestrial radiation.
+    See :py:func:`~xclim.indices._helpers.extraterrestrial_solar_radiation`.
 
     References
     ----------
@@ -1944,7 +1944,6 @@ def mean_radiant_temperature(
         of the solar zenith angle is calculated. If "sunlit", the cosine of the
         solar zenith angle is calculated during the sunlit period of each interval.
         If "instant", the instantaneous cosine of the solar zenith angle is calculated.
-        This is necessary if mrt is not None.
 
     Returns
     -------
@@ -1969,46 +1968,27 @@ def mean_radiant_temperature(
     rlus = convert_units_to(rlus, "W m-2")
 
     dates = rsds.time
-    hours = ((dates - dates.dt.floor("D")).dt.seconds / 3600).assign_attrs(units="h")
-
     lat = _gather_lat(rsds)
     lon = _gather_lon(rsds)
+    dec = solar_declination(dates)
 
-    decimal_year = datetime_to_decimal_year(times=dates, calendar=dates.dt.calendar)
-    day_angle = ((decimal_year % 1) * 2 * np.pi).assign_attrs(units="rad")
-    dec = solar_declination(day_angle)
     if stat == "sunlit":
-        interval = (dates.diff("time").dt.seconds / 3600).reindex(
-            time=dates.time, method="bfill"
-        )
         csza_i = cosine_of_solar_zenith_angle(
-            declination=dec,
-            lat=lat,
-            lon=lon,
-            hours=hours,
-            interval=interval,
-            stat="interval",
+            dates, dec, lat, lon=lon, stat="average", sunlit=False
         )
         csza_s = cosine_of_solar_zenith_angle(
-            declination=dec, lat=lat, lon=lon, hours=hours, interval=interval, stat=stat
+            dates, dec, lat, lon=lon, stat="average", sunlit=True
         )
     elif stat == "instant":
-        tc = time_correction_for_solar_angle(day_angle)
+        tc = time_correction_for_solar_angle(dates)
         csza = cosine_of_solar_zenith_angle(
-            declination=dec,
-            lat=lat,
-            lon=lon,
-            time_correction=tc,
-            hours=hours,
-            stat="instant",
+            dates, dec, lat, lon=lon, time_correction=tc, stat="instant"
         )
         csza_i = csza.copy()
         csza_s = csza.copy()
     elif stat == "average":
         csza = cosine_of_solar_zenith_angle(
-            declination=dec,
-            lat=lat,
-            stat="average",
+            dates, dec, lat, stat="average", sunlit=False
         )
         csza_i = csza.copy()
         csza_s = csza.copy()