Repository: Ouranosinc/xclim

Issue:
`subset_shape` with large input data
* xclim version: 0.15.0
* Python version: 3.8
* Operating System: ubuntu

### Description

Calling the `subset_shape` function on a large input dataset crashes when the dataset doesn't fit in memory. Not with a large polygon subset, but a large source dataset. So I want to subset a small portion of a large file.


### What I Did

I was using xclim through finch, the source data is the daily bccaqv2 datasets. 

### What I Received

```
finch         | INFO:PYWPS:Subsetting file 1 of 48 
finch         | Traceback (most recent call last):
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/pywps/app/Process.py", line 243, in _run_process
finch         |     self.handler(wps_request, wps_response)  # the user must update the wps_response.
finch         |   File "/code/finch/processes/wps_base.py", line 36, in _handler_wrapper
finch         |     return self.wrapped_handler(request, response)
finch         |   File "/code/finch/processes/wps_ensemble_indices_polygon.py", line 74, in _handler
finch         |     return ensemble_common_handler(self, request, response, finch_subset_shape)
finch         |   File "/code/finch/processes/ensemble_utils.py", line 528, in ensemble_common_handler
finch         |     subsetted_files = subset_function(
finch         |   File "/code/finch/processes/subset.py", line 234, in finch_subset_shape
finch         |     process_threaded(_subset, netcdf_inputs)
finch         |   File "/code/finch/processes/utils.py", line 274, in process_threaded
finch         |     outputs = [function(r) for r in inputs]
finch         |   File "/code/finch/processes/utils.py", line 274, in <listcomp>
finch         |     outputs = [function(r) for r in inputs]
finch         |   File "/code/finch/processes/subset.py", line 219, in _subset
finch         |     subsetted = subset_shape(
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/xclim/subset.py", line 208, in func_checker
finch         |     final = func(*formatted_args, **kwargs)
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/xclim/subset.py", line 465, in subset_shape
finch         |     ds_copy = copy.deepcopy(ds)
finch         |   File "/opt/conda/envs/finch/lib/python3.8/copy.py", line 153, in deepcopy
finch         |     y = copier(memo)
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/xarray/core/dataset.py", line 1181, in __deepcopy__
finch         |     return self.copy(deep=True)
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/xarray/core/dataset.py", line 1073, in copy
finch         |     variables = {k: v.copy(deep=deep) for k, v in self._variables.items()}
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/xarray/core/dataset.py", line 1073, in <dictcomp>
finch         |     variables = {k: v.copy(deep=deep) for k, v in self._variables.items()}
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/xarray/core/variable.py", line 885, in copy
finch         |     data = np.array(data)
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/xarray/core/indexing.py", line 677, in __array__
finch         |     self._ensure_cached()
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/xarray/core/indexing.py", line 674, in _ensure_cached
finch         |     self.array = NumpyIndexingAdapter(np.asarray(self.array))
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/numpy/core/_asarray.py", line 85, in asarray
finch         |     return array(a, dtype, copy=False, order=order)
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/xarray/core/indexing.py", line 653, in __array__
finch         |     return np.asarray(self.array, dtype=dtype)
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/numpy/core/_asarray.py", line 85, in asarray
finch         |     return array(a, dtype, copy=False, order=order)
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/xarray/core/indexing.py", line 557, in __array__
finch         |     return np.asarray(array[self.key], dtype=None)
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/numpy/core/_asarray.py", line 85, in asarray
finch         |     return array(a, dtype, copy=False, order=order)
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/xarray/coding/variables.py", line 72, in __array__
finch         |     return self.func(self.array)
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/xarray/coding/variables.py", line 218, in _scale_offset_decoding
finch         |     data = np.array(data, dtype=dtype, copy=True)
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/xarray/coding/variables.py", line 72, in __array__
finch         |     return self.func(self.array)
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/xarray/coding/variables.py", line 138, in _apply_mask
finch         |     data = np.asarray(data, dtype=dtype)
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/numpy/core/_asarray.py", line 85, in asarray
finch         |     return array(a, dtype, copy=False, order=order)
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/xarray/core/indexing.py", line 557, in __array__
finch         |     return np.asarray(array[self.key], dtype=None)
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/xarray/backends/netCDF4_.py", line 72, in __getitem__
finch         |     return indexing.explicit_indexing_adapter(
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/xarray/core/indexing.py", line 837, in explicit_indexing_adapter
finch         |     result = raw_indexing_method(raw_key.tuple)
finch         |   File "/opt/conda/envs/finch/lib/python3.8/site-packages/xarray/backends/netCDF4_.py", line 85, in _getitem
finch         |     array = getitem(original_array, key)
finch         |   File "netCDF4/_netCDF4.pyx", line 4389, in netCDF4._netCDF4.Variable.__getitem__
finch         | MemoryError: Unable to allocate 55.9 GiB for an array with shape (1068, 510, 55115) and data type int16
```


Discussion:

Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index 1660b080..d6263236 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -8,6 +8,8 @@ History
 * CRS definitions of projected DataSets are now written to file according to Climate and Forecast-convention standards
 * Add utilities to merge attributes and update history in xclim.core.formatting
 * Ensembles : Allow alignment of datasets with same frequency but different offsets
+* Bug fixes in run_length for run-with-dates methods when the date is not found in the run.
+* Remove deepcopy from subset.subset_shape to improve memory usage.
 
 0.15.x (2020-03-12)
 -------------------
diff --git a/xclim/indices/_multivariate.py b/xclim/indices/_multivariate.py
index b8014fe0..f861fd08 100644
--- a/xclim/indices/_multivariate.py
+++ b/xclim/indices/_multivariate.py
@@ -107,7 +107,7 @@ def cold_spell_duration_index(
 
     below = tasmin < thresh
 
-    return below.resample(time=freq).apply(
+    return below.resample(time=freq).map(
         rl.windowed_run_count, window=window, dim="time"
     )
 
@@ -565,7 +565,7 @@ def heat_wave_frequency(
 
     cond = (tasmin > thresh_tasmin) & (tasmax > thresh_tasmax)
     group = cond.resample(time=freq)
-    return group.apply(rl.windowed_run_events, window=window, dim="time")
+    return group.map(rl.windowed_run_events, window=window, dim="time")
 
 
 @declare_units(
@@ -636,7 +636,7 @@ def heat_wave_max_length(
 
     cond = (tasmin > thresh_tasmin) & (tasmax > thresh_tasmax)
     group = cond.resample(time=freq)
-    max_l = group.apply(rl.longest_run, dim="time")
+    max_l = group.map(rl.longest_run, dim="time")
     return max_l.where(max_l >= window, 0)
 
 
@@ -691,7 +691,7 @@ def heat_wave_total_length(
 
     cond = (tasmin > thresh_tasmin) & (tasmax > thresh_tasmax)
     group = cond.resample(time=freq)
-    return group.apply(rl.windowed_run_count, args=(window,), dim="time")
+    return group.map(rl.windowed_run_count, args=(window,), dim="time")
 
 
 @declare_units("", pr="[precipitation]", prsn="[precipitation]", tas="[temperature]")
@@ -1374,7 +1374,7 @@ def warm_spell_duration_index(
 
     above = tasmax > thresh
 
-    return above.resample(time=freq).apply(
+    return above.resample(time=freq).map(
         rl.windowed_run_count, window=window, dim="time"
     )
 
diff --git a/xclim/indices/_simple.py b/xclim/indices/_simple.py
index 210e4e42..a0cce9fb 100644
--- a/xclim/indices/_simple.py
+++ b/xclim/indices/_simple.py
@@ -414,7 +414,7 @@ def consecutive_frost_days(tasmin: xarray.DataArray, freq: str = "AS-JUL"):
     if fu != tu:
         frz = units.convert(frz, fu, tu)
     group = (tasmin < frz).resample(time=freq)
-    return group.apply(rl.longest_run, dim="time")
+    return group.map(rl.longest_run, dim="time")
 
 
 @declare_units("days", tasmin="[temperature]")
diff --git a/xclim/indices/_threshold.py b/xclim/indices/_threshold.py
index 2e3b1be6..68e25bb0 100644
--- a/xclim/indices/_threshold.py
+++ b/xclim/indices/_threshold.py
@@ -89,7 +89,7 @@ def cold_spell_days(
     over = tas < t
     group = over.resample(time=freq)
 
-    return group.apply(rl.windowed_run_count, window=window, dim="time")
+    return group.map(rl.windowed_run_count, window=window, dim="time")
 
 
 @declare_units("mm/day", pr="[precipitation]", thresh="[precipitation]")
@@ -223,7 +223,7 @@ def maximum_consecutive_wet_days(
     thresh = convert_units_to(thresh, pr, "hydro")
 
     group = (pr > thresh).resample(time=freq)
-    return group.apply(rl.longest_run, dim="time")
+    return group.map(rl.longest_run, dim="time")
 
 
 @declare_units("C days", tas="[temperature]", thresh="[temperature]")
@@ -458,7 +458,7 @@ def growing_season_length(
     )
 
 
-@declare_units("days", tas="[temperature]", thresh="[temperature]")
+@declare_units("", tas="[temperature]", thresh="[temperature]")
 def last_spring_frost(
     tas: xarray.DataArray,
     thresh: str = "0 degC",
@@ -533,7 +533,7 @@ def heat_wave_index(
     over = tasmax > thresh
     group = over.resample(time=freq)
 
-    return group.apply(rl.windowed_run_count, window=window, dim="time")
+    return group.map(rl.windowed_run_count, window=window, dim="time")
 
 
 @declare_units("C days", tas="[temperature]", thresh="[temperature]")
@@ -628,7 +628,7 @@ def hot_spell_max_length(
 
     cond = tasmax > thresh_tasmax
     group = cond.resample(time=freq)
-    max_l = group.apply(rl.longest_run, dim="time")
+    max_l = group.map(rl.longest_run, dim="time")
     return max_l.where(max_l >= window, 0)
 
 
@@ -686,7 +686,7 @@ def hot_spell_frequency(
 
     cond = tasmax > thresh_tasmax
     group = cond.resample(time=freq)
-    return group.apply(rl.windowed_run_events, window=window, dim="time")
+    return group.map(rl.windowed_run_events, window=window, dim="time")
 
 
 @declare_units("days", tasmin="[temperature]", thresh="[temperature]")
@@ -901,7 +901,7 @@ def maximum_consecutive_dry_days(
     t = convert_units_to(thresh, pr, "hydro")
     group = (pr < t).resample(time=freq)
 
-    return group.apply(rl.longest_run, dim="time")
+    return group.map(rl.longest_run, dim="time")
 
 
 @declare_units("days", tasmin="[temperature]", thresh="[temperature]")
@@ -944,7 +944,7 @@ def maximum_consecutive_frost_free_days(
     t = convert_units_to(thresh, tasmin)
     group = (tasmin > t).resample(time=freq)
 
-    return group.apply(rl.longest_run, dim="time")
+    return group.map(rl.longest_run, dim="time")
 
 
 @declare_units("days", tasmax="[temperature]", thresh="[temperature]")
@@ -987,7 +987,7 @@ def maximum_consecutive_tx_days(
     t = convert_units_to(thresh, tasmax)
     group = (tasmax > t).resample(time=freq)
 
-    return group.apply(rl.longest_run, dim="time")
+    return group.map(rl.longest_run, dim="time")
 
 
 @declare_units("[area]", sic="[]", area="[area]", thresh="[]")
diff --git a/xclim/indices/generic.py b/xclim/indices/generic.py
index 5ec2020c..d16a9173 100644
--- a/xclim/indices/generic.py
+++ b/xclim/indices/generic.py
@@ -69,7 +69,7 @@ def select_resample_op(da: xr.DataArray, op, freq: str = "YS", **indexer):
     if isinstance(op, str):
         return getattr(r, op)(dim="time", keep_attrs=True)
 
-    return r.apply(op)
+    return r.map(op)
 
 
 def doymax(da: xr.DataArray):
diff --git a/xclim/indices/run_length.py b/xclim/indices/run_length.py
index e5c9b55e..3e0866da 100644
--- a/xclim/indices/run_length.py
+++ b/xclim/indices/run_length.py
@@ -248,7 +248,10 @@ def first_run(
         if isinstance(coord, str):
             crd = getattr(crd.dt, coord)
 
-        return lazy_indexing(crd, out)
+        out = lazy_indexing(crd, out)
+
+    if dim in out:
+        out = out.drop_vars(dim)
 
     return out
 
@@ -321,30 +324,20 @@ def run_length_with_date(
     """
     include_date = datetime.strptime(date, "%m-%d").timetuple().tm_yday
 
-    def _crld(rlda, rlwindow, rldoy, rldim):
-        mid_index = np.where(rlda.time.dt.dayofyear == rldoy)[0]
-        if (
-            mid_index.size == 0
-        ):  # The date is not within the group. Happens at boundaries.
-            all_nans = xr.full_like(rlda.isel(time=0), np.nan)
-            all_nans.attrs = {}
-            return all_nans
-        end = first_run(
-            (~rlda).where(rlda.time >= rlda.time[mid_index][0]),
-            window=rlwindow,
-            dim=rldim,
-        )
-        beg = first_run(rlda, window=rlwindow, dim=rldim)
-        sl = end - beg
-        sl = xr.where(
-            beg.isnull() & end.notnull(), 0, sl
-        )  # If series is never triggered
-        sl = xr.where(
-            beg.notnull() & end.isnull(), rlda.time.size - beg, sl
-        )  # If series is not ended by end of resample time frequency
-        return sl.where(sl >= 0)
-
-    return _crld(rlda=da, rlwindow=window, rldoy=include_date, rldim=dim)
+    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.
+        return xr.full_like(da.isel(time=0), np.nan, float).drop_vars("time")
+
+    end = first_run(
+        (~da).where(da.time >= da.time[mid_index][0]), window=window, dim=dim,
+    )
+    beg = first_run(da, window=window, dim=dim)
+    sl = end - beg
+    sl = xr.where(beg.isnull() & end.notnull(), 0, sl)  # If series is never triggered
+    sl = xr.where(
+        beg.notnull() & end.isnull(), da.time.size - beg, sl
+    )  # If series is not ended by end of resample time frequency
+    return sl.where(sl >= 0)
 
 
 def run_end_after_date(
@@ -378,29 +371,21 @@ def run_end_after_date(
     """
     after_date = datetime.strptime(date, "%m-%d").timetuple().tm_yday
 
-    def _cred(rlda, rlwindow, rldoy, rldim, rlcoord):
-        mid_idx = np.where(rlda.time.dt.dayofyear == rldoy)[0]
-        if (
-            mid_idx.size == 0
-        ):  # The date is not within the group. Happens at boundaries.
-            all_nans = xr.full_like(rlda.isel(time=0), np.nan)
-            all_nans.attrs = {}
-            return all_nans
-        end = first_run(
-            (~rlda).where(rlda.time >= rlda.time[mid_idx][0]),
-            window=rlwindow,
-            dim=rldim,
-            coord=rlcoord,
-        )
-        beg = first_run(
-            rlda.where(rlda.time < rlda.time[mid_idx][0]), window=rlwindow, dim=rldim,
-        )
-        end = xr.where(
-            end.isnull() & beg.notnull(), rlda.time.isel(time=-1).dt.dayofyear, end
-        )
-        return end.where(beg.notnull())
-
-    return _cred(rlda=da, rlwindow=window, rldoy=after_date, rldim=dim, rlcoord=coord)
+    mid_idx = np.where(da.time.dt.dayofyear == after_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")
+
+    end = first_run(
+        (~da).where(da.time >= da.time[mid_idx][0]),
+        window=window,
+        dim=dim,
+        coord=coord,
+    )
+    beg = first_run(da.where(da.time < da.time[mid_idx][0]), window=window, dim=dim)
+    end = xr.where(
+        end.isnull() & beg.notnull(), da.time.isel(time=-1).dt.dayofyear, end
+    )
+    return end.where(beg.notnull())
 
 
 def last_run_before_date(
@@ -434,20 +419,12 @@ def last_run_before_date(
     """
     before_date = datetime.strptime(date, "%m-%d").timetuple().tm_yday
 
-    def _clrbd(rlda, rlwindow, rlbefore_date, rldim, rlcoord):
-        mid_idx = np.where(rlda.time.dt.dayofyear == rlbefore_date)[0]
-        if (
-            mid_idx.size == 0
-        ):  # The date is not within the group. Happens at boundaries.
-            all_nans = xr.full_like(rlda.isel(time=0), np.nan)
-            all_nans.attrs = dict()
-            return all_nans
-        run = rlda.where(rlda.time <= rlda.time[mid_idx][0])
-        return last_run(run, window=rlwindow, dim=rldim, coord=rlcoord)
-
-    return _clrbd(
-        rlda=da, rlwindow=window, rlbefore_date=before_date, rldim=dim, rlcoord=coord
-    )
+    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")
+
+    run = da.where(da.time <= da.time[mid_idx][0])
+    return last_run(run, window=window, dim=dim, coord=coord)
 
 
 def rle_1d(
diff --git a/xclim/subset.py b/xclim/subset.py
index 778876c5..c83e4a31 100644
--- a/xclim/subset.py
+++ b/xclim/subset.py
@@ -1,4 +1,3 @@
-import copy
 import logging
 import warnings
 from functools import wraps
@@ -458,11 +457,10 @@ def subset_shape(
     """
     wgs84 = CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
 
-    # TODO : edge case using polygon splitting decorator touches original ds when subsetting?
     if isinstance(ds, xarray.DataArray):
-        ds_copy = copy.deepcopy(ds._to_temp_dataset())
+        ds_copy = ds._to_temp_dataset()
     else:
-        ds_copy = copy.deepcopy(ds)
+        ds_copy = ds.copy()
 
     if isinstance(shape, gpd.GeoDataFrame):
         poly = shape.copy()


Test patch:
diff --git a/tests/test_run_length.py b/tests/test_run_length.py
index 072f5052..53342644 100644
--- a/tests/test_run_length.py
+++ b/tests/test_run_length.py
@@ -37,15 +37,15 @@ class TestLongestRun:
         )
         values[1:11] = 1
         da = xr.DataArray(values != 0, coords={"time": time}, dims="time")
-        lt = da.resample(time="M").apply(rl.longest_run_ufunc)
+        lt = da.resample(time="M").map(rl.longest_run_ufunc)
         assert lt[0] == 10
         np.testing.assert_array_equal(lt[1:], 0)
 
         # n-dim version versus ufunc
         da3d = xr.open_dataset(self.nc_pr).pr[:, 40:50, 50:68] != 0
-        lt_orig = da3d.resample(time="M").apply(rl.longest_run_ufunc)
+        lt_orig = da3d.resample(time="M").map(rl.longest_run_ufunc)
         # override 'auto' usage of ufunc for small number of gridpoints
-        lt_Ndim = da3d.resample(time="M").apply(
+        lt_Ndim = da3d.resample(time="M").map(
             rl.longest_run, dim="time", ufunc_1dim=False
         )
         np.testing.assert_array_equal(lt_orig, lt_Ndim)
@@ -57,7 +57,7 @@ class TestLongestRun:
         )
         values[0:10] = 1
         da = xr.DataArray(values != 0, coords={"time": time}, dims="time")
-        lt = da.resample(time="M").apply(rl.longest_run_ufunc)
+        lt = da.resample(time="M").map(rl.longest_run_ufunc)
         assert lt[0] == 10
         np.testing.assert_array_equal(lt[1:], 0)
 
@@ -65,9 +65,9 @@ class TestLongestRun:
         da3d = xr.open_dataset(self.nc_pr).pr[:, 40:50, 50:68] * 0
         da3d[0:10,] = da3d[0:10,] + 1
         da3d = da3d == 1
-        lt_orig = da3d.resample(time="M").apply(rl.longest_run_ufunc)
+        lt_orig = da3d.resample(time="M").map(rl.longest_run_ufunc)
         # override 'auto' usage of ufunc for small number of gridpoints
-        lt_Ndim = da3d.resample(time="M").apply(
+        lt_Ndim = da3d.resample(time="M").map(
             rl.longest_run, dim="time", ufunc_1dim=False
         )  # override 'auto' for small
         np.testing.assert_array_equal(lt_orig, lt_Ndim)
@@ -80,7 +80,7 @@ class TestLongestRun:
         values[-10:] = 1
         da = xr.DataArray(values != 0, coords={"time": time}, dims="time")
 
-        lt = da.resample(time="M").apply(rl.longest_run_ufunc)
+        lt = da.resample(time="M").map(rl.longest_run_ufunc)
         assert lt[-1] == 10
         np.testing.assert_array_equal(lt[:-1], 0)
 
@@ -88,8 +88,8 @@ class TestLongestRun:
         da3d = xr.open_dataset(self.nc_pr).pr[:, 40:50, 50:68] * 0
         da3d[-10:,] = da3d[-10:,] + 1
         da3d = da3d == 1
-        lt_orig = da3d.resample(time="M").apply(rl.longest_run_ufunc)
-        lt_Ndim = da3d.resample(time="M").apply(
+        lt_orig = da3d.resample(time="M").map(rl.longest_run_ufunc)
+        lt_Ndim = da3d.resample(time="M").map(
             rl.longest_run, dim="time", ufunc_1dim=False
         )
         np.testing.assert_array_equal(lt_orig, lt_Ndim)
@@ -101,14 +101,14 @@ class TestLongestRun:
         )
         da = xr.DataArray(values != 0, coords={"time": time}, dims="time")
 
-        lt = da.resample(time="M").apply(rl.longest_run_ufunc)
+        lt = da.resample(time="M").map(rl.longest_run_ufunc)
         np.testing.assert_array_equal(lt, da.resample(time="M").count(dim="time"))
 
         # n-dim version versus ufunc
         da3d = xr.open_dataset(self.nc_pr).pr[:, 40:50, 50:68] * 0 + 1
         da3d = da3d == 1
-        lt_orig = da3d.resample(time="M").apply(rl.longest_run_ufunc)
-        lt_Ndim = da3d.resample(time="M").apply(
+        lt_orig = da3d.resample(time="M").map(rl.longest_run_ufunc)
+        lt_Ndim = da3d.resample(time="M").map(
             rl.longest_run, dim="time", ufunc_1dim=False
         )
         np.testing.assert_array_equal(lt_orig, lt_Ndim)
@@ -121,7 +121,7 @@ class TestLongestRun:
         )
         da = xr.DataArray(values != 0, coords={"time": time}, dims="time")
 
-        lt = da.resample(time="M").apply(rl.longest_run_ufunc)
+        lt = da.resample(time="M").map(rl.longest_run_ufunc)
         n = da.resample(time="M").count(dim="time")
         np.testing.assert_array_equal(lt[0], n[0])
         np.testing.assert_array_equal(lt[1], 26)
@@ -130,8 +130,8 @@ class TestLongestRun:
         da3d = xr.open_dataset(self.nc_pr).pr[:, 40:50, 50:68] * 0 + 1
         da3d[35,] = da3d[35,] + 1
         da3d = da3d == 1
-        lt_orig = da3d.resample(time="M").apply(rl.longest_run_ufunc)
-        lt_Ndim = da3d.resample(time="M").apply(
+        lt_orig = da3d.resample(time="M").map(rl.longest_run_ufunc)
+        lt_Ndim = da3d.resample(time="M").map(
             rl.longest_run, dim="time", ufunc_1dim=False
         )
         np.testing.assert_array_equal(lt_orig, lt_Ndim)
@@ -186,8 +186,8 @@ class TestWindowedRunEvents:
 
         # n-dim version versus ufunc
         da3d = xr.open_dataset(self.nc_pr).pr[:, 40:50, 50:68] != 0
-        lt_orig = da3d.resample(time="M").apply(rl.windowed_run_events_ufunc, window=4)
-        lt_Ndim = da3d.resample(time="M").apply(
+        lt_orig = da3d.resample(time="M").map(rl.windowed_run_events_ufunc, window=4)
+        lt_Ndim = da3d.resample(time="M").map(
             rl.windowed_run_events, window=4, dim="time", ufunc_1dim=False
         )
         np.testing.assert_array_equal(lt_orig, lt_Ndim)
@@ -204,8 +204,8 @@ class TestWindowedRunCount:
 
         # n-dim version versus ufunc
         da3d = xr.open_dataset(self.nc_pr).pr[:, 40:50, 50:68] != 0
-        lt_orig = da3d.resample(time="M").apply(rl.windowed_run_count_ufunc, window=4)
-        lt_Ndim = da3d.resample(time="M").apply(
+        lt_orig = da3d.resample(time="M").map(rl.windowed_run_count_ufunc, window=4)
+        lt_Ndim = da3d.resample(time="M").map(
             rl.windowed_run_count, window=4, dim="time", ufunc_1dim=False
         )
         np.testing.assert_array_equal(lt_orig, lt_Ndim)
@@ -305,3 +305,20 @@ class TestRunsWithDates:
             runs, window=1, date=date, dim="time", coord=coord
         )
         np.testing.assert_array_equal(np.mean(out.load()), expected)
+
+    @pytest.mark.parametrize(
+        "func",
+        [rl.last_run_before_date, rl.run_end_after_date, rl.run_length_with_date],
+    )
+    @pytest.mark.parametrize("use_dask", [True, False])
+    def test_run_with_dates_no_date(self, tas_series, use_dask, func):
+        t = np.ones(90)
+        tas = tas_series(t, start="2000-01-01")
+        runs = xr.concat((tas, tas), dim="dim0")
+        runs = runs == 1
+
+        if use_dask:
+            runs = runs.chunk({"time": 10, "dim0": 1})
+
+        out = runs.resample(time="MS").map(func, window=1, date="07-01", dim="time")
+        assert out.isnull().all()