Repository: Ouranosinc/xclim

Issue:
Run calculations with dates are not calendar-aware
* xclim version: 0.20
* Python version: 3.8
* Operating System: Linux

### Description
<!--Describe what you were trying to get done.
Tell us what happened, what went wrong, and what you expected to happen.-->
Functions of `xclim.indices.run_length` that accept "date" arguments all transform those dates into day-of-year integers, but do so without considering the data's calendar. 

### What I Did
<!--Paste the command(s) you ran and the output.
If there was a crash, please include the traceback below.-->
Synthetic example with 250-days run of warm temp starting in 03-01.
```python3
for cal in ['standard', '365_day', '366_day']:    
    time = xr.cftime_range('2004-01-01', end='2005-12-31', freq='D', calendar=cal)
    
    tas = np.zeros(time.size)
    start = np.where((time.day == 1)  & (time.month == 3))[0]
    # default cal 31 + 29, 
    tas[start[0]:start[0] + 250] = 5
    tas[start[1]:start[1] + 250] = 5
    tas = xr.DataArray(tas, coords={'time': time}, dims=('time',))
    print(cal, (tas > 0).resample(time='AS-MAR').map(rl.first_run_after_date, date='03-01', window=2).values)
```
outputs:
```
standard [nan nan 60.]  # Expected:  [nan, 61, 60]
365_day [nan 60. 60.]  # Expected: [nan, 60, 60]
366_day [nan nan nan]  # Expected [nan, 61, 61]
```
What happens is that `03-01` is translated using the default calendar for a non-leap year : `doy = 60`. When resampling with `time='AS-MAR'`, non leap years have 60 as first doy but leap years have 61. Thus, the "after_date" condition masks the whole sequence and the method returns nan.

Discussion:

Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index 67715d9a..04b7dcd9 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -18,6 +18,7 @@ New features and enhancements
 Bug fixes
 ~~~~~~~~~
 * Fixed a bug in the attributes of `frost_season_length`.
+* `indices.run_length` methods using dates now respect the array's calendar.
 
 Internal changes
 ~~~~~~~~~~~~~~~~~
diff --git a/xclim/indices/run_length.py b/xclim/indices/run_length.py
index a65e1191..237a8e35 100644
--- a/xclim/indices/run_length.py
+++ b/xclim/indices/run_length.py
@@ -349,14 +349,11 @@ def season_length(
     not_da = (~da).where(da.time.copy(data=np.arange(da.time.size)) >= beg.fillna(0))
     if date is not None:
         # Mask also values after "date"
-        include_date = datetime.strptime(date, "%m-%d").timetuple().tm_yday
-
-        mid_index = np.where(da.time.dt.dayofyear == include_date)[0]
-        if (
-            mid_index.size == 0
-        ):  # The date is not within the group. Happens at boundaries.
+        mid_idx = index_of_date(da.time, date, max_idxs=1)
+        if mid_idx.size == 0:
+            # The date is not within the group. Happens at boundaries.
             return xr.full_like(da.isel(time=0), np.nan, float).drop_vars("time")
-        not_da = not_da.where(da.time >= da.time[mid_index][0])
+        not_da = not_da.where(da.time >= da.time[mid_idx][0])
 
     end = first_run(
         not_da,
@@ -369,7 +366,7 @@ def season_length(
         beg.notnull() & end.isnull(), da.time.size - beg, sl
     )  # If series is not ended by end of resample time frequency
     if date is not None:
-        sl = sl.where(beg < mid_index.squeeze())
+        sl = sl.where(beg < mid_idx.squeeze())
     sl = xr.where(beg.isnull() & end.notnull(), 0, sl)  # If series is never triggered
     return sl
 
@@ -405,9 +402,7 @@ def run_end_after_date(
     xr.DataArray
       Index (or coordinate if `coord` is not False) of last item in last valid run. Returns np.nan if there are no valid run.
     """
-    after_date = datetime.strptime(date, "%m-%d").timetuple().tm_yday
-
-    mid_idx = np.where(da.time.dt.dayofyear == after_date)[0]
+    mid_idx = index_of_date(da.time, date, max_idxs=1)
     if mid_idx.size == 0:  # The date is not within the group. Happens at boundaries.
         return xr.full_like(da.isel(time=0), np.nan, float).drop_vars("time")
 
@@ -421,7 +416,7 @@ def run_end_after_date(
     end = xr.where(
         end.isnull() & beg.notnull(), da.time.isel(time=-1).dt.dayofyear, end
     )
-    return end.where(beg.notnull())
+    return end.where(beg.notnull()).drop_vars("time")
 
 
 def first_run_after_date(
@@ -453,9 +448,7 @@ def first_run_after_date(
     out : xr.DataArray
       Index (or coordinate if `coord` is not False) of first item in the first valid run. Returns np.nan if there are no valid run.
     """
-    after_date = datetime.strptime(date, "%m-%d").timetuple().tm_yday
-
-    mid_idx = np.where(da.time.dt.dayofyear == after_date)[0]
+    mid_idx = index_of_date(da.time, date, max_idxs=1)
     if mid_idx.size == 0:  # The date is not within the group. Happens at boundaries.
         return xr.full_like(da.isel(time=0), np.nan, float).drop_vars("time")
 
@@ -496,9 +489,8 @@ def last_run_before_date(
     xr.DataArray
       Index (or coordinate if `coord` is not False) of last item in last valid run. Returns np.nan if there are no valid run.
     """
-    before_date = datetime.strptime(date, "%m-%d").timetuple().tm_yday
+    mid_idx = index_of_date(da.time, date)
 
-    mid_idx = np.where(da.time.dt.dayofyear == before_date)[0]
     if mid_idx.size == 0:  # The date is not within the group. Happens at boundaries.
         return xr.full_like(da.isel(time=0), np.nan, float).drop_vars("time")
 
@@ -797,3 +789,35 @@ def lazy_indexing(da: xr.DataArray, index: xr.DataArray, dim=None):
         dask="parallelized",
         output_dtypes=[da.dtype],
     )
+
+
+def index_of_date(time: xr.DataArray, date: str, max_idxs: Optional[int] = None):
+    """Get the index of a date in a time array.
+
+    Parameters
+    ----------
+    time : xr.DataArray
+      An array of datetime values, any calendar.
+    date : str
+      A string in the "yyyy-mm-dd" or "mm-dd" format.
+
+    Returns
+    -------
+    ndarray
+      1D array of integers, indexes of `date` in `time`.
+    """
+    try:
+        date = datetime.strptime(date, "%Y-%m-%d")
+        year_cond = time.dt.year == date.year
+    except ValueError:
+        date = datetime.strptime(date, "%m-%d")
+        year_cond = True
+
+    idxs = np.where(
+        year_cond & (time.dt.month == date.month) & (time.dt.day == date.day)
+    )[0]
+    if max_idxs is not None and idxs.size > max_idxs:
+        raise ValueError(
+            f"More than {max_idxs} instance of date {date} found in the coordinate array."
+        )
+    return idxs


Test patch:
diff --git a/tests/test_indices.py b/tests/test_indices.py
index 185843ba..85c10884 100644
--- a/tests/test_indices.py
+++ b/tests/test_indices.py
@@ -404,7 +404,7 @@ class TestGrowingSeasonEnd:
             ("2000-07-10", "2001-01-01", "07-01", np.nan),  # End happens before start
             ("2000-06-15", "2000-07-15", "07-01", 198),  # Normal case
             ("2000-06-15", "2000-07-25", "07-15", 208),  # PCC Case
-            ("2000-06-15", "2000-07-15", "10-01", 274),  # Late mid_date
+            ("2000-06-15", "2000-07-15", "10-01", 275),  # Late mid_date
             ("2000-06-15", "2000-07-15", "01-10", np.nan),  # Early mid_date
         ],
     )
diff --git a/tests/test_run_length.py b/tests/test_run_length.py
index 6543d5d3..35b6263c 100644
--- a/tests/test_run_length.py
+++ b/tests/test_run_length.py
@@ -255,7 +255,7 @@ class TestRunsWithDates:
             ("07-01", 210, 70),
             ("07-01", 190, 50),
             ("04-01", 150, np.NaN),  # date falls early
-            ("11-01", 150, 164),  # date ends late
+            ("11-01", 150, 165),  # date ends late
             (None, 150, 10),  # no date, real length
         ],
     )
@@ -284,7 +284,7 @@ class TestRunsWithDates:
             ("dayofyear", "07-01", 210, 211),
             (False, "07-01", 190, 190),
             ("dayofyear", "04-01", 150, np.NaN),  # date falls early
-            ("dayofyear", "11-01", 150, 305),  # date ends late
+            ("dayofyear", "11-01", 150, 306),  # date ends late
         ],
     )
     @pytest.mark.parametrize("use_dask", [True, False])
@@ -307,7 +307,7 @@ class TestRunsWithDates:
             ("dayofyear", "07-01", 210, 211),
             (False, "07-01", 190, 190),
             ("dayofyear", "04-01", False, np.NaN),  # no run
-            ("dayofyear", "11-01", 150, 305),  # run already started
+            ("dayofyear", "11-01", 150, 306),  # run already started
         ],
     )
     @pytest.mark.parametrize("use_dask", [True, False])
@@ -332,8 +332,8 @@ class TestRunsWithDates:
     @pytest.mark.parametrize(
         "coord,date,end,expected",
         [
-            ("dayofyear", "07-01", 210, 182),
-            (False, "07-01", 190, 181),
+            ("dayofyear", "07-01", 210, 183),
+            (False, "07-01", 190, 182),
             ("dayofyear", "04-01", 150, np.NaN),  # date falls early
             ("dayofyear", "11-01", 150, 150),  # date ends late
         ],
@@ -373,6 +373,55 @@ class TestRunsWithDates:
         out = runs.resample(time="MS").map(func, window=1, date="07-01", dim="time")
         assert out.isnull().all()
 
+    @pytest.mark.parametrize(
+        "calendar,expected",
+        [("standard", [61, 60]), ("365_day", [60, 60]), ("366_day", [61, 61])],
+    )
+    def test_run_with_dates_different_calendars(self, calendar, expected):
+        time = xr.cftime_range(
+            "2004-01-01", end="2005-12-31", freq="D", calendar=calendar
+        )
+        tas = np.zeros(time.size)
+        start = np.where((time.day == 1) & (time.month == 3))[0]
+        tas[start[0] : start[0] + 250] = 5
+        tas[start[1] : start[1] + 250] = 5
+        tas = xr.DataArray(tas, coords={"time": time}, dims=("time",))
+        out = (
+            (tas > 0)
+            .resample(time="AS-MAR")
+            .map(rl.first_run_after_date, date="03-01", window=2)
+        )
+        np.testing.assert_array_equal(out.values[1:], expected)
+
+        out = (
+            (tas > 0)
+            .resample(time="AS-MAR")
+            .map(rl.season_length, date="03-02", window=2)
+        )
+        np.testing.assert_array_equal(out.values[1:], [250, 250])
+
+        out = (
+            (tas > 0)
+            .resample(time="AS-MAR")
+            .map(rl.run_end_after_date, date="03-03", window=2)
+        )
+        np.testing.assert_array_equal(out.values[1:], np.array(expected) + 250)
+
+        out = (
+            (tas > 0)
+            .resample(time="AS-MAR")
+            .map(rl.last_run_before_date, date="03-02", window=2)
+        )
+        np.testing.assert_array_equal(out.values[1:], np.array(expected) + 1)
+
+    @pytest.mark.parametrize(
+        "func", [rl.first_run_after_date, rl.season_length, rl.run_end_after_date]
+    )
+    def test_too_many_dates(self, func, tas_series):
+        tas = tas_series(np.zeros(730), start="2000-01-01")
+        with pytest.raises(ValueError, match="More than 1 instance of date"):
+            func((tas == 0), date="03-01", window=5)
+
 
 @pytest.mark.parametrize("use_dask", [True, False])
 def test_lazy_indexing_special_cases(use_dask):