Repository: Ouranosinc/xclim

Issue:
frequency_analysis crashes if there are empty grid cells
* xclim version: 0.11.2
* Python version: 3.6
* Operating System:

### Description
I have a Dataset on a rotated grid (EC-EARTH_HIRHAM5, from CORDEX-NA). I've made a subset for my region, and on the borders of that subset, some grid cells are empty (meaning that all values are NaN).

Following a call to frequency_analysis, the code crashes if I try to save the data to NetCDF, or even with a simple ".values" call. I receive the following error:
ValueError: zero-size array to reduction operation minimum which has no identity


### What I Did
```
import xarray as xr
from xclim import subset
from xclim import generic as xcgen

ds = xr.open_mfdataset("/path/to/rotated/pole/data")  <--- EC-EARTH_HIRHAM5 hourly precipitation for the 1981-2010 period
ds_subset = subset.subset_bbox(ds, lat_bnds=lat_bnds, lon_bnds=lon_bnds)
da_freq = xcgen.frequency_analysis(ds_subset["pr"], mode="max", t=2, dist="genextreme", window=6, freq="YS")
da_freq.isel(return_period=0).values  <-- crashes
```

If I use the following code, everything works fine (except that I get data in locations that should be NaN):
```
import xarray as xr
from xclim import subset
from xclim import generic as xcgen

ds = xr.open_mfdataset("/path/to/rotated/pole/data")
ds_subset = subset.subset_bbox(ds, lat_bnds=lat_bnds, lon_bnds=lon_bnds)
ds_subset["pr"] = ds_subset["pr"].fillna(value=0)  <-- this is the new line
da_freq = xcgen.frequency_analysis(ds_subset["pr"], mode="max", t=2, dist="genextreme", window=6, freq="YS")
da_freq.isel(return_period=0).values  
```


Discussion:

Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index ced91a46..a42ac318 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -15,6 +15,7 @@ History
 * Reformatted packaging configuration to pure Py3 wheel that ignore tests and test data
 * Now officially supporting Python3.8
 * Enhancement to precip_accumulation() to allow estimated amounts solid (or liquid) phase precipitation
+* Bugfix for frequency analysis choking on time series with NaNs only.
 
 0.10.x-beta (2019-06-18)
 ------------------------
diff --git a/xclim/generic.py b/xclim/generic.py
index 0fee5a9f..14d6868a 100644
--- a/xclim/generic.py
+++ b/xclim/generic.py
@@ -102,31 +102,42 @@ def fit(da, dist="norm"):
     """
     # Get the distribution
     dc = get_dist(dist)
+    shape_params = [] if dc.shapes is None else dc.shapes.split(",")
+    dist_params = shape_params + ["loc", "scale"]
 
-    # Fit the parameters (lazy computation)
-    data = dask.array.apply_along_axis(
-        dc.fit, da.get_axis_num("time"), da.dropna("time", how="all")
-    )
+    # Fit the parameters.
+    # This would also be the place to impose constraints on the series minimum length if needed.
+    def fitfunc(arr):
+        """Fit distribution parameters."""
+        x = np.ma.masked_invalid(arr).compressed()
 
-    # Count the number of values used for the fit.
-    # n = arr.count(dim='time')
+        # Return NaNs if array is empty.
+        if len(x) <= 1:
+            return [np.nan,] * len(dist_params)
 
-    # Create a view to a DataArray with the desired dimensions to copy them over to the parameter array.
-    mean = da.mean(dim="time", keep_attrs=True)
+        # Fill with NaNs if one of the parameters is NaN
+        params = dc.fit(x)
+        if np.isnan(params).any():
+            params[:] = np.nan
 
-    # Create coordinate for the distribution parameters
-    coords = dict(mean.coords.items())
-    coords["dparams"] = ([] if dc.shapes is None else dc.shapes.split(",")) + [
-        "loc",
-        "scale",
-    ]
+        return params
 
-    # TODO: add time and time_bnds coordinates (Low will work on this)
-    # time.attrs['climatology'] = 'climatology_bounds'
-    # coords['time'] =
-    # coords['climatology_bounds'] =
+    # xarray.apply_ufunc does not yet support multiple outputs with dask parallelism.
+    data = dask.array.apply_along_axis(fitfunc, da.get_axis_num("time"), da)
 
-    out = xr.DataArray(data=data, coords=coords, dims=("dparams",) + mean.dims)
+    # Count the number of values used for the fit.
+    # n = da.notnull().count(dim='time')
+
+    # Coordinates for the distribution parameters
+    coords = dict(da.coords.items())
+    coords.pop("time")
+    coords["dparams"] = dist_params
+
+    # Dimensions for the distribution parameters
+    dims = ["dparams",] + list(da.dims)
+    dims.remove("time")
+
+    out = xr.DataArray(data=data, coords=coords, dims=dims)
     out.attrs = da.attrs
     out.attrs["original_name"] = getattr(da, "standard_name", "")
     out.attrs[
@@ -260,7 +271,7 @@ def frequency_analysis(da, mode, t, dist, window=1, freq=None, **indexer):
     freq = freq or default_freq(**indexer)
 
     # Extract the time series of min or max over the period
-    sel = select_resample_op(da, op=mode, freq=freq, **indexer).dropna(dim="time")
+    sel = select_resample_op(da, op=mode, freq=freq, **indexer)
 
     # Frequency analysis
     return fa(sel, t, dist, mode)


Test patch:
diff --git a/tests/test_generic.py b/tests/test_generic.py
index 38b040b1..67b0a572 100644
--- a/tests/test_generic.py
+++ b/tests/test_generic.py
@@ -42,6 +42,36 @@ class TestFA(object):
         q0 = lognorm.ppf(1 - 1.0 / T, *p0)
         np.testing.assert_array_equal(q[0, 0, 0], q0)
 
+    def test_fit_nan(self):
+        da = self.da.copy()
+        da[0, 0, 0] = np.nan
+        out_nan = generic.fit(da, "lognorm")
+        out_censor = generic.fit(da[1:], "lognorm")
+        np.testing.assert_array_equal(
+            out_nan.values[:, 0, 0], out_censor.values[:, 0, 0]
+        )
+
+    def test_empty(self):
+        da = self.da.copy()
+        da[:, 0, 0] = np.nan
+        out = generic.fit(da, "lognorm").values
+        assert np.isnan(out[:, 0, 0]).all()
+
+
+class TestFrequencyAnalysis:
+    def test_simple(self, ndq_series):
+        q = ndq_series.copy()
+        q[:, 0, 0] = np.nan
+        out = generic.frequency_analysis(
+            q, mode="max", t=2, dist="genextreme", window=6, freq="YS"
+        )
+        assert out.dims == ("return_period", "x", "y")
+        assert out.shape == (1, 2, 3)
+        v = out.values
+        assert v.shape == (1, 2, 3)
+        assert np.isnan(v[:, 0, 0])
+        assert ~np.isnan(v[:, 1, 1])
+
 
 class TestSelectResampleOp:
     def test_month(self, q_series):
diff --git a/tests/test_streamflow.py b/tests/test_streamflow.py
index 41906f4a..b1aa1397 100644
--- a/tests/test_streamflow.py
+++ b/tests/test_streamflow.py
@@ -28,6 +28,14 @@ class Test_FA:
         )
         assert out.shape == (1, 2, 3)
 
+    def test_empty(self, ndq_series):
+        q = ndq_series.copy()
+        q[:, 0, 0] = np.nan
+        out = streamflow.freq_analysis(
+            q, mode="max", t=2, dist="genextreme", window=6, freq="YS"
+        )
+        assert np.isnan(out.values[:, 0, 0]).all()
+
 
 class TestStats:
     def test_simple(self, ndq_series):