Remove faulty default in MRT and UTCI (#1501)

<!--Please ensure the PR fulfills the following requirements! -->
<!-- If this is your first PR, make sure to add your details to the
AUTHORS.rst! -->
### Pull Request Checklist:
- [x] This PR addresses an already opened issue (for bug fixes /
features)
    - This PR fixes #1496
- [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?
I believe #1399 changed the meanings of argument `stat` which made the
default option of `mean_radiant_temperature` incorrect. This PR removes
the "average" option and changes the default to 'sunlit'.

### Does this PR introduce a breaking change?
Yes, outputs will change. However, I think previous output was
incorrect. An announcement has been added to the CHANGES.

### Other information:

.github/workflows/main.yml

@@ -61,7 +61,7 @@ jobs:
     strategy:
       matrix:
         include:
-          - tox-env: "py39"
+          - tox-env: "py39-coverage"
             python-version: "3.9"
     steps:
       - uses: actions/[email protected]

CHANGES.rst

@@ -10,6 +10,10 @@ New indicators
 ^^^^^^^^^^^^^^
 * ``xclim.indices.snw_storm_days`` computes the number of days with snowfall amount accumulation above a threshold. (:pull:`1505`).
 
+Announcements
+^^^^^^^^^^^^^
+* The default mechanism for computing the Mean Radiant Temperature, a part of the Universal Thermal Climate Index (UTCI) was broken in xclim 0.44 and 0.45. This has been fixed by changing the default.
+
 New features and enhancements
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 * Add ``wind_power_potential`` to estimate the potential for wind power production given wind speed at the turbine hub height and turbine specifications, along with  ``wind_profile`` to estimate the wind speed at different heights based on wind speed at a reference height. (:issue:`1458`, :pull:`1471`)
@@ -36,6 +40,7 @@ Bug fixes
 * Indicators ``snd_season_length`` and ``snw_season_length`` will return 0 instead of NaN if all inputs have a (non-NaN) zero snow depth (or water-equivalent thickness). (:pull:`1492`, :issue:`1491`)
 * Fixed a bug in the `pytest` configuration that could prevent testing data caching from occurring in systems where the platform-dependent cache directory is not found in the user's home. (:issue:`1468`, :pull:`1473`).
 * Fix ``xclim.core.dataflags.data_flags`` variable name generation (:pull:`1507`).
+* Remove nonsensical `stat='average'` option for ``mean_radiant_temperature``. (:issue:`1496`, :pull:`1501`).
 
 Breaking changes
 ^^^^^^^^^^^^^^^^

tests/test_atmos.py

@@ -587,7 +587,7 @@ def test_universal_thermal_climate_index(self, atmosds):
             rsus=rsus,
             rlds=rlds,
             rlus=rlus,
-            stat="average",
+            stat="sunlit",
         )
 
         np.testing.assert_allclose(utci.isel(time=0), utci_exp, rtol=1e-03)
@@ -605,15 +605,13 @@ def test_mean_radiant_temperature(self, atmosds):
         # Expected values
         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")
         mrt_ins = atmos.mean_radiant_temperature(rsds, rsus, rlds, rlus, stat="instant")
-        mrt_avg = atmos.mean_radiant_temperature(rsds, rsus, rlds, rlus, stat="average")
+
         rtol = 1e-03
         np.testing.assert_allclose(mrt_sun.isel(time=0), exp_sun, rtol=rtol)
         np.testing.assert_allclose(mrt_ins.isel(time=0), exp_ins, rtol=rtol)
-        np.testing.assert_allclose(mrt_avg.isel(time=0), exp_avg, rtol=rtol)
 
 
 class TestLateFrostDays:

tests/test_indices.py

@@ -3423,9 +3423,7 @@ def test_universal_thermal_climate_index(
     np.testing.assert_allclose(utci, utci_exp, rtol=1e-03)
 
 
-@pytest.mark.parametrize(
-    "stat,expected", [("sunlit", 295.0), ("instant", 294.9), ("average", 295.1)]
-)
+@pytest.mark.parametrize("stat,expected", [("sunlit", 295.0), ("instant", 294.9)])
 def test_mean_radiant_temperature(
     rsds_series,
     rsus_series,

xclim/indices/_conversion.py

@@ -1775,7 +1775,7 @@ def universal_thermal_climate_index(
     rsus: xr.DataArray = None,
     rlds: xr.DataArray = None,
     rlus: xr.DataArray = None,
-    stat: str = "average",
+    stat: str = "sunlit",
     mask_invalid: bool = True,
 ) -> xr.DataArray:
     r"""Universal thermal climate index (UTCI).
@@ -1805,12 +1805,10 @@ def universal_thermal_climate_index(
     rlus : xr.DataArray, optional
         Surface Upwelling Longwave Radiation
         This is necessary if mrt is not None.
-    stat : {'average', 'instant', 'sunlit'}
-        Which statistic to apply. If "average", the average of the cosine of the
-        solar zenith angle is calculated. If "instant", the instantaneous cosine
+    stat : {'instant', 'sunlit'}
+        Which statistic to apply. If "instant", the instantaneous cosine
         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.
     mask_invalid: bool
         If True (default), UTCI values are NaN where any of the inputs are outside
@@ -1874,8 +1872,7 @@ def universal_thermal_climate_index(
 
 def _fdir_ratio(
     dates: xr.DataArray,
-    csza_i: xr.DataArray,
-    csza_s: xr.DataArray,
+    csza: xr.DataArray,
     rsds: xr.DataArray,
 ) -> xr.DataArray:
     r"""Return ratio of direct solar radiation.
@@ -1887,10 +1884,8 @@ def _fdir_ratio(
     ----------
     dates : xr.DataArray
         Series of dates and time of day
-    csza_i : xr.DataArray
-        Cosine of the solar zenith angle during each interval
-    csza_s : xr.DataArray
-        Cosine of the solar zenith angle during the sunlit period of each interval
+    csza : xr.DataArray
+        Cosine of the solar zenith angle during the sunlit period of each interval or at an instant
     rsds : xr.DataArray
         Surface Downwelling Shortwave Radiation
 
@@ -1908,12 +1903,12 @@ def _fdir_ratio(
     :cite:cts:`liljegren_modeling_2008,kong_explicit_2022`
     """
     d = distance_from_sun(dates)
-    s_star = rsds * ((1367 * csza_s * (d ** (-2))) ** (-1))
+    s_star = rsds * ((1367 * csza * (d ** (-2))) ** (-1))
     s_star = xr.where(s_star > 0.85, 0.85, s_star)
     fdir_ratio = np.exp(3 - 1.34 * s_star - 1.65 * (s_star ** (-1)))
     fdir_ratio = xr.where(fdir_ratio > 0.9, 0.9, fdir_ratio)
     return xr.where(
-        (fdir_ratio <= 0) | (csza_i <= np.cos(89.5 / 180 * np.pi)) | (rsds <= 0),
+        (fdir_ratio <= 0) | (csza <= np.cos(89.5 / 180 * np.pi)) | (rsds <= 0),
         0,
         fdir_ratio,
     )
@@ -1927,7 +1922,7 @@ def mean_radiant_temperature(
     rsus: xr.DataArray,
     rlds: xr.DataArray,
     rlus: xr.DataArray,
-    stat: str = "average",
+    stat: str = "sunlit",
 ) -> xr.DataArray:
     r"""Mean radiant temperature.
 
@@ -1943,12 +1938,10 @@ def mean_radiant_temperature(
         Surface Downwelling Longwave Radiation
     rlus : xr.DataArray
         Surface Upwelling Longwave Radiation
-    stat : {'average', 'instant', 'sunlit'}
-        Which statistic to apply. If "average", the average of the cosine of the
-        solar zenith angle is calculated. If "instant", the instantaneous cosine
+    stat : {'instant', 'sunlit'}
+        Which statistic to apply. If "instant", the instantaneous cosine
         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.
 
     Returns
     -------
@@ -1978,38 +1971,27 @@ def mean_radiant_temperature(
     dec = solar_declination(dates)
 
     if stat == "sunlit":
-        csza_i = cosine_of_solar_zenith_angle(
-            dates, dec, lat, lon=lon, stat="average", sunlit=False
-        )
-        csza_s = cosine_of_solar_zenith_angle(
+        csza = cosine_of_solar_zenith_angle(
             dates, dec, lat, lon=lon, stat="average", sunlit=True
         )
     elif stat == "instant":
         tc = time_correction_for_solar_angle(dates)
         csza = cosine_of_solar_zenith_angle(
             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(
-            dates, dec, lat, stat="average", sunlit=False
-        )
-        csza_i = csza.copy()
-        csza_s = csza.copy()
     else:
         raise NotImplementedError(
-            "Argument 'stat' must be one of 'average', 'instant' or 'sunlit'."
+            "Argument 'stat' must be one of 'instant' or 'sunlit'."
         )
 
-    fdir_ratio = _fdir_ratio(dates, csza_i, csza_s, rsds)
+    fdir_ratio = _fdir_ratio(dates, csza, rsds)
 
     rsds_direct = fdir_ratio * rsds
     rsds_diffuse = rsds - rsds_direct
 
-    gamma = np.arcsin(csza_i)
+    gamma = np.arcsin(csza)
     fp = 0.308 * np.cos(gamma * 0.988 - (gamma**2 / 50000))
-    i_star = xr.where(csza_s > 0.001, rsds_direct / csza_s, 0)
+    i_star = xr.where(csza > 0.001, rsds_direct / csza, 0)
 
     mrt = np.power(
         (

xclim/indices/helpers.py

@@ -197,7 +197,7 @@ def cosine_of_solar_zenith_angle(
     lat: Quantified,
     lon: Quantified = "0 °",
     time_correction: xr.DataArray = None,
-    stat: str = "integral",
+    stat: str = "average",
     sunlit: bool = False,
 ) -> xr.DataArray:
     """Cosine of the solar zenith angle.
@@ -222,10 +222,10 @@ def cosine_of_solar_zenith_angle(
     time_correction : xr.DataArray, optional
         Time correction for solar angle. See :py:func:`time_correction_for_solar_angle`
         This is necessary if stat is "instant".
-    stat : {'integral', 'average', 'instant'}
+    stat : {'average', 'integral', 'instant'}
         Which daily statistic to return.
-        If "integral", this returns the integral of the cosine of the zenith angle
         If "average", this returns the average of the cosine of the zenith angle
+        If "integral", this returns the integral of the cosine of the zenith angle
         If "instant", this returns the instantaneous cosine of the zenith angle
     sunlit: bool
         If True, only the sunlit part of the interval is considered in the integral or average.