Snow season length is 0 if there is 0 snow - fix solar_zenith_angle t…

…ypo (#1492)

<!--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 #1491 and fixes #1493
- [x] Tests for the changes have been added (for bug fixes / features)
- [ ] (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?
`snd_season_length` and `snw_season_length` got implemented with the
same check as their `*_season_start` and `*_season_end` counter parts,
with a mask where : `not at_least_n_valid(snw.where(snw > 0), n=1,
freq=freq)`.

However, this means the output is NaN if all inputs are 0. This makes
sense in the DOY case (can't a have a start/end date if there's no
season), but it doesn't in the "length" case. Instead, I think (as does
the issue raiser) an all-0 snow timeseries simply means a season length
of 0. Thus, I removed the `.where(snw > 0)` part of the test for those
two indicators.

EDIT: I also slipped in another fix for
`xc.indices.helpers.cosine_of_solar_zenith_angle`. There was a typo in
the `xr.apply_ufunc` call.

### Does this PR introduce a breaking change?
Yes because it will change the output of two indicators, but I think
this is for the best. The previous behaviour was not documented as such
in the doc, I believe the new one is actually more intuitive. No need
for a warning, IMO.

.github/workflows/main.yml

@@ -75,6 +75,8 @@ jobs:
         run: tox -e ${{ matrix.tox-env }}
         env:
           GITHUB_TOKEN: ${{ secrets.GITHUB_TOKEN }}
+          COVERALLS_FLAG_NAME: run-{{ matrix.tox-env }}
+          COVERALLS_SERVICE_NAME: github
 
   test-pypi:
     needs: lint

CHANGES.rst

@@ -4,7 +4,7 @@ Changelog
 
 v0.46.0 (unreleased)
 --------------------
-Contributors to this version: Éric Dupuis (:user:`coxipi`), Trevor James Smith (:user:`Zeitsperre`), David Huard (:user:`huard`).
+Contributors to this version: Éric Dupuis (:user:`coxipi`), Trevor James Smith (:user:`Zeitsperre`), David Huard (:user:`huard`) and Pascal Bourgault (:user:`aulemahal`).
 
 New features and enhancements
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
@@ -21,6 +21,7 @@ Bug fixes
 * Calling a ``sdba.map_blocks``-wrapped function with data chunked along the reduced dimensions will raise an error. This forbids chunking the trained dataset along the distribution dimensions, for example. (:issue:`1481`, :pull:`1482`).
 * Optimization of indicators ``huglin_index`` and ``biologically_effective_degree_days`` when used with dask and flox. As a side effect, the indice functions (i.e. under ``xc.indices``) no longer mask incomplete periods. The indicators' output is unchanged under the default "check_missing" setting (:issue:`1494`, :pull:`1495`).
 * Fixed ``xclim.indices.run_length.lazy_indexing`` which would sometimes trigger the loading of auxiliary coordinates. (:issue:`1483`, :pull:`1484`).
+* 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`).
 
 Breaking changes

tests/test_indices.py

@@ -2682,20 +2682,21 @@ def test_nan_slices(self, snd_series, snw_series):
 
 
 class TestSnowCover:
-    def test_snow_season_length(self, snd_series, snw_series):
-        a = np.ones(366) / 100.0
-        a[10:20] = 0.3
+    @pytest.mark.parametrize("length", [0, 10])
+    def test_snow_season_length(self, snd_series, snw_series, length):
+        a = np.zeros(366)
+        a[10 : 10 + length] = 0.3
         snd = snd_series(a)
         # kg m-2 = 1000 kg m-3 * 1 m
         snw = snw_series(1000 * a)
 
         out = xci.snd_season_length(snd)
         assert len(out) == 2
-        assert out[0] == 10
+        assert out[0] == length
 
         out = xci.snw_season_length(snw)
         assert len(out) == 2
-        assert out[0] == 10
+        assert out[0] == length
 
     def test_continous_snow_season_start(self, snd_series, snw_series):
         a = np.arange(366) / 100.0

tests/test_snow.py

@@ -25,11 +25,14 @@ def test_simple(self, snd_series):
 
 
 class TestSnowWaterCoverDuration:
-    def test_simple(self, snw_series):
-        snw = snw_series(np.ones(110) * 1000, start="2001-01-01")
+    @pytest.mark.parametrize(
+        "factor,exp", ([1000, [31, 28, 31, np.NaN]], [0, [0, 0, 0, np.NaN]])
+    )
+    def test_simple(self, snw_series, factor, exp):
+        snw = snw_series(np.ones(110) * factor, start="2001-01-01")
         out = land.snw_season_length(snw, freq="M")
         assert out.units == "days"
-        np.testing.assert_array_equal(out, [31, 28, 31, np.nan])
+        np.testing.assert_array_equal(out, exp)
 
 
 class TestContinuousSnowDepthCoverStartEnd:

xclim/indices/_threshold.py

@@ -2069,7 +2069,7 @@ def snd_season_length(
     xarray.DataArray, [time]
         Number of days where snow depth is greater than or equal to threshold.
     """
-    valid = at_least_n_valid(snd.where(snd > 0), n=1, freq=freq)
+    valid = at_least_n_valid(snd, n=1, freq=freq)
     thresh = convert_units_to(thresh, snd)
     out = threshold_count(snd, op, thresh, freq)
     return to_agg_units(out, snd, "count").where(~valid)
@@ -2106,7 +2106,7 @@ def snw_season_length(
     xarray.DataArray, [time]
         Number of days where snow water is greater than or equal to threshold.
     """
-    valid = at_least_n_valid(snw.where(snw > 0), n=1, freq=freq)
+    valid = at_least_n_valid(snw, n=1, freq=freq)
     thresh = convert_units_to(thresh, snw)
     out = threshold_count(snw, op, thresh, freq)
     return to_agg_units(out, snw, "count").where(~valid)

xclim/indices/helpers.py

@@ -296,7 +296,7 @@ def cosine_of_solar_zenith_angle(
         _wrap_radians(h_e),
         stat == "average",
         input_core_dims=[[]] * 6,
-        dask="parallel",
+        dask="parallelized",
     )