Repository: Ouranosinc/xclim

Issue:
Add fit indicator
### Description

The `fit` function has no related indicator instance.



Discussion:

Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index 13f89f90..710cdc1a 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -6,6 +6,7 @@ History
 --------------
 * Migrated to a `major.minor.patch-release` semantic versioning system.
 * Removed attributes in netCDF output from Indicators that are not in the CF-convention.
+* Added `fit` indicator to fit the parameters of a distribution to a series.
 
 0.10-beta (06-06-2019)
 ----------------------
diff --git a/setup.cfg b/setup.cfg
index 161509ba..3aede051 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,15 +1,15 @@
 [bumpversion]
-current_version = 0.10.0-beta
+current_version = 0.10.1-beta
 commit = True
 tag = True
 parse = (?P<major>\d+)\.(?P<minor>\d+).(?P<patch>\d+)(\-(?P<release>[a-z]+))?
-serialize =
+serialize = 
 	{major}.{minor}.{patch}-{release}
 	{major}.{minor}.{patch}
 
 [bumpversion:part:release]
 optional_value = gamma
-values =
+values = 
 	beta
 	gamma
 
@@ -25,14 +25,14 @@ replace = __version__ = '{new_version}'
 universal = 1
 
 [flake8]
-exclude =
+exclude = 
 	.git,
 	docs,
 	build,
 	.eggs,
 max-line-length = 120
 max-complexity = 10
-ignore =
+ignore = 
 	F401
 	F403
 
@@ -42,6 +42,6 @@ test = pytest
 [tool:pytest]
 collect_ignore = ['setup.py']
 addopts = --verbose
-filterwarnings =
+filterwarnings = 
 	ignore::UserWarning
 
diff --git a/setup.py b/setup.py
index 5d55cb57..ce599b4d 100644
--- a/setup.py
+++ b/setup.py
@@ -64,6 +64,6 @@ setup(
     tests_require=test_requirements,
     extras_require={'docs': docs_requirements},
     url='https://github.com/Ouranosinc/xclim',
-    version='0.10.0-beta',
+    version='0.10.1-beta',
     zip_safe=False,
 )
diff --git a/xclim/__init__.py b/xclim/__init__.py
index 22588fb6..36f66bfe 100644
--- a/xclim/__init__.py
+++ b/xclim/__init__.py
@@ -11,7 +11,7 @@ import sys
 
 __author__ = """Travis Logan"""
 __email__ = 'logan.travis@ouranos.ca'
-__version__ = '0.10.0-beta'
+__version__ = '0.10.1-beta'
 
 
 def build_module(name, objs, doc='', source=None, mode='ignore'):
diff --git a/xclim/generic.py b/xclim/generic.py
index fa422ac2..1aee24f8 100644
--- a/xclim/generic.py
+++ b/xclim/generic.py
@@ -80,12 +80,12 @@ def doymin(da):
     return out
 
 
-def fit(arr, dist='norm'):
+def fit(da, dist='norm'):
     """Fit an array to a univariate distribution along the time dimension.
 
     Parameters
     ----------
-    arr : xarray.DataArray
+    da : xarray.DataArray
       Time series to be fitted along the time dimension.
     dist : str
       Name of the univariate distribution, such as beta, expon, genextreme, gamma, gumbel_r, lognorm, norm
@@ -95,18 +95,23 @@ def fit(arr, dist='norm'):
     -------
     xarray.DataArray
       An array of distribution parameters fitted using the method of Maximum Likelihood.
+
+    Notes
+    -----
+    Coordinates for which all values are NaNs will be dropped before fitting the distribution. If the array
+    still contains NaNs, the distribution parameters will be returned as NaNs.
     """
     # Get the distribution
     dc = get_dist(dist)
 
     # Fit the parameters (lazy computation)
-    data = dask.array.apply_along_axis(dc.fit, arr.get_axis_num('time'), arr)
+    data = dask.array.apply_along_axis(dc.fit, da.get_axis_num('time'), da.dropna('time', how='all'))
 
     # Count the number of values used for the fit.
     # n = arr.count(dim='time')
 
     # Create a view to a DataArray with the desired dimensions to copy them over to the parameter array.
-    mean = arr.mean(dim='time', keep_attrs=True)
+    mean = da.mean(dim='time', keep_attrs=True)
 
     # Create coordinate for the distribution parameters
     coords = dict(mean.coords.items())
@@ -118,25 +123,23 @@ def fit(arr, dist='norm'):
     # coords['climatology_bounds'] =
 
     out = xr.DataArray(data=data, coords=coords, dims=(u'dparams',) + mean.dims)
-    out.attrs = arr.attrs
-    out.attrs['original_name'] = getattr(arr, 'standard_name', '')
-    out.attrs['standard_name'] = '{0} distribution parameters'.format(dist)
-    out.attrs['long_name'] = '{0} distribution parameters for {1}'.format(dist, getattr(arr, 'standard_name', ''))
+    out.attrs = da.attrs
+    out.attrs['original_name'] = getattr(da, 'standard_name', '')
+    out.attrs['description'] = \
+        'Parameters of the {0} distribution fitted over {1}'.format(dist, getattr(da, 'standard_name', ''))
     out.attrs['estimator'] = 'Maximum likelihood'
-    out.attrs['cell_methods'] = (out.attrs.get('cell_methods', '') + ' time: fit').strip()
+    out.attrs['scipy_dist'] = dist
     out.attrs['units'] = ''
-    msg = '\nData fitted with {0} statistical distribution using a Maximum Likelihood Estimator'
-    out.attrs['history'] = out.attrs.get('history', '') + msg.format(dist)
-
+    # out.name = 'params'
     return out
 
 
-def fa(arr, t, dist='norm', mode='high'):
+def fa(da, t, dist='norm', mode='high'):
     """Return the value corresponding to the given return period.
 
     Parameters
     ----------
-    arr : xarray.DataArray
+    da : xarray.DataArray
       Maximized/minimized input data with a `time` dimension.
     t : int or sequence
       Return period. The period depends on the resolution of the input data. If the input array's resolution is
@@ -158,7 +161,7 @@ def fa(arr, t, dist='norm', mode='high'):
     dc = get_dist(dist)
 
     # Fit the parameters of the distribution
-    p = fit(arr, dist)
+    p = fit(da, dist)
 
     # Create a lambda function to facilitate passing arguments to dask. There is probably a better way to do this.
     if mode in ['max', 'high']:
@@ -190,9 +193,9 @@ def fa(arr, t, dist='norm', mode='high'):
     out = xr.DataArray(data=data, coords=coords, dims=dims)
     out.attrs = p.attrs
     out.attrs['standard_name'] = '{0} quantiles'.format(dist)
-    out.attrs['long_name'] = '{0} return period values for {1}'.format(dist, getattr(arr, 'standard_name', ''))
+    out.attrs['long_name'] = '{0} return period values for {1}'.format(dist, getattr(da, 'standard_name', ''))
     out.attrs['cell_methods'] = (out.attrs.get('cell_methods', '') + ' dparams: ppf').strip()
-    out.attrs['units'] = arr.attrs.get('units', '')
+    out.attrs['units'] = da.attrs.get('units', '')
     out.attrs['mode'] = mode
     out.attrs['history'] = out.attrs.get('history', '') + "Compute values corresponding to return periods."
 
diff --git a/xclim/streamflow.py b/xclim/streamflow.py
index 29353098..c4e8d05a 100644
--- a/xclim/streamflow.py
+++ b/xclim/streamflow.py
@@ -34,13 +34,19 @@ class Stats(Streamflow):
         return reduce(np.logical_or, miss)
 
 
-# We need to disable the missing value check because the output here is not a time series.
+# Disable the missing value check because the output here is not a time series.
 class FA(Streamflow):
     def missing(self, *args, **kwds):
         """Return whether an output is considered missing or not."""
         return False
 
 
+# Disable the daily checks because the inputs are period extremas.
+class Fit(FA):
+    def validate(self, da):
+        pass
+
+
 base_flow_index = Streamflow(identifier='base_flow_index',
                              units='',
                              long_name="Base flow index",
@@ -62,6 +68,16 @@ stats = Stats(identifier='stats',
               compute=generic.select_resample_op)
 
 
+fit = Fit(identifier='fit',
+          var_name='params',
+          units='',
+          standard_name='{dist} parameters',
+          long_name="{dist} distribution parameters",
+          description="Parameters of the {dist} distribution fitted ",
+          cell_methods='time: fit',
+          compute=generic.fit)
+
+
 doy_qmax = Streamflow(identifier='doy_qmax',
                       var_name='q{indexer}_doy_qmax',
                       long_name='Day of the year of the maximum over {indexer}',
diff --git a/xclim/utils.py b/xclim/utils.py
index 6f6ec015..9a167322 100644
--- a/xclim/utils.py
+++ b/xclim/utils.py
@@ -648,7 +648,7 @@ class Indicator():
 
         # Fill in missing metadata from the doc
         meta = parse_doc(self.compute.__doc__)
-        for key in ['abstract', 'title', 'long_name', 'notes', 'references']:
+        for key in ['abstract', 'title', 'notes', 'references']:
             setattr(self, key, getattr(self, key) or meta.get(key, ''))
 
     def __call__(self, *args, **kwds):
@@ -687,7 +687,7 @@ class Indicator():
 
         attrs['history'] += '[{:%Y-%m-%d %H:%M:%S}] {}: {}{}'.format(
             dt.datetime.now(), vname, self.identifier, ba.signature.replace(parameters=cp.values()))
-        attrs['cell_methods'] += out_attrs.pop('cell_methods')
+        attrs['cell_methods'] += out_attrs.pop('cell_methods', '')
         attrs.update(out_attrs)
 
         # Assume the first arguments are always the DataArray.
@@ -721,7 +721,7 @@ class Indicator():
         """CF-Convention attributes of the output value."""
         names = ['standard_name', 'long_name', 'units', 'cell_methods', 'description', 'comment',
                  'references']
-        return {k: getattr(self, k, '') for k in names}
+        return {k: getattr(self, k) for k in names if getattr(self, k)}
 
     def json(self, args=None):
         """Return a dictionary representation of the class.


Test patch:
diff --git a/tests/test_generic.py b/tests/test_generic.py
index d3591d64..ea4ca001 100644
--- a/tests/test_generic.py
+++ b/tests/test_generic.py
@@ -32,6 +32,7 @@ class TestFA(object):
         # Check that we can reuse the parameters with scipy distributions
         cdf = lognorm.cdf(.99, *p.values)
         assert cdf.shape == (self.nx, self.ny)
+        assert p.attrs['estimator'] == 'Maximum likelihood'
 
     def test_fa(self):
         T = 10
diff --git a/tests/test_streamflow.py b/tests/test_streamflow.py
index e8afd450..90b101ab 100644
--- a/tests/test_streamflow.py
+++ b/tests/test_streamflow.py
@@ -40,6 +40,14 @@ class TestStats():
         np.testing.assert_array_equal(out[2].isnull(), True)
 
 
+class TestFit():
+
+    def test_simple(self, ndq_series):
+        ts = streamflow.stats(ndq_series, freq='YS', op='max')
+        p = streamflow.fit(ts, dist='gumbel_r')
+        assert p.attrs['estimator'] == 'Maximum likelihood'
+
+
 def test_qdoy_max(ndq_series, q_series):
     out = streamflow.doy_qmax(ndq_series, freq='YS', season='JJA')
     assert out.attrs['units'] == ''