diff --git a/CHANGELOG.rst b/CHANGELOG.rst
index 9b2d1484..e4e726a3 100644
--- a/CHANGELOG.rst
+++ b/CHANGELOG.rst
@@ -12,6 +12,7 @@ Features
   observations are probability distributions ``p`` or cumulative
   distributionss ``c``. See :py:func:`~xskillscore.rps` docstrings and doctests for
   examples. (:pr:`300`) `Aaron Spring`_
+- Added slope of linear fit :py:func:`~xskillscore.linslope`. `Ray Bell`_
 
 Internal Changes
 ~~~~~~~~~~~~~~~~
diff --git a/ci/docs_notebooks.yml b/ci/docs_notebooks.yml
index 540ad2bf..6e4b78db 100644
--- a/ci/docs_notebooks.yml
+++ b/ci/docs_notebooks.yml
@@ -25,6 +25,6 @@ dependencies:
   - sphinx_rtd_theme
   - pip
   - pip:
-    - sphinx_autosummary_accessors
-    # Install latest version of xskillscore.
-    - -e ..
+      - sphinx_autosummary_accessors
+      # Install latest version of xskillscore.
+      - -e ..
diff --git a/ci/minimum-tests.yml b/ci/minimum-tests.yml
index ce70b61b..69f47170 100644
--- a/ci/minimum-tests.yml
+++ b/ci/minimum-tests.yml
@@ -23,4 +23,4 @@ dependencies:
   - pytest-xdist
   - pip
   - pip:
-    - -e ..
+      - -e ..
diff --git a/docs/source/api.rst b/docs/source/api.rst
index a59cbc43..a351298a 100644
--- a/docs/source/api.rst
+++ b/docs/source/api.rst
@@ -24,6 +24,7 @@ Correlation Metrics
     spearman_r_eff_p_value
     effective_sample_size
     r2
+    linslope
 
 Distance Metrics
 ~~~~~~~~~~~~~~~~
diff --git a/docs/source/api/xskillscore.linslope.rst b/docs/source/api/xskillscore.linslope.rst
new file mode 100644
index 00000000..f714a499
--- /dev/null
+++ b/docs/source/api/xskillscore.linslope.rst
@@ -0,0 +1,6 @@
+xskillscore.linslope
+====================
+
+.. currentmodule:: xskillscore
+
+.. autofunction:: linslope
diff --git a/docs/source/quick-start.ipynb b/docs/source/quick-start.ipynb
index d31375e7..d0e10e6e 100644
--- a/docs/source/quick-start.ipynb
+++ b/docs/source/quick-start.ipynb
@@ -66,6 +66,7 @@
     "* Spearman Correlation p value (`spearman_r_p_value`)\n",
     "* Spearman Correlation effective p value (`spearman_r_eff_p_value`)\n",
     "* Effective Sample Size (`effective_sample_size`)\n",
+    "* Slope of Linear Fit (`linslope`)\n",
     "* Coefficient of Determination (`r2`)\n",
     "\n",
     "### Distance-Based\n",
diff --git a/xskillscore/__init__.py b/xskillscore/__init__.py
index ff9018c5..9366a13c 100644
--- a/xskillscore/__init__.py
+++ b/xskillscore/__init__.py
@@ -7,6 +7,7 @@
 from .core.contingency import Contingency
 from .core.deterministic import (
     effective_sample_size,
+    linslope,
     mae,
     mape,
     me,
diff --git a/xskillscore/core/accessor.py b/xskillscore/core/accessor.py
index 9bfc41be..31a70f5b 100644
--- a/xskillscore/core/accessor.py
+++ b/xskillscore/core/accessor.py
@@ -2,6 +2,7 @@
 
 from .deterministic import (
     effective_sample_size,
+    linslope,
     mae,
     mape,
     me,
@@ -46,6 +47,11 @@ def _in_ds(self, x):
         else:
             return self._obj[x]
 
+    def linslope(self, a, b, *args, **kwargs):
+        a = self._in_ds(a)
+        b = self._in_ds(b)
+        return linslope(a, b, *args, **kwargs)
+
     def pearson_r(self, a, b, *args, **kwargs):
         a = self._in_ds(a)
         b = self._in_ds(b)
diff --git a/xskillscore/core/deterministic.py b/xskillscore/core/deterministic.py
index 201f8662..9bae56d3 100644
--- a/xskillscore/core/deterministic.py
+++ b/xskillscore/core/deterministic.py
@@ -4,6 +4,7 @@
 
 from .np_deterministic import (
     _effective_sample_size,
+    _linslope,
     _mae,
     _mape,
     _me,
@@ -27,21 +28,22 @@
 )
 
 __all__ = [
+    "effective_sample_size",
+    "linslope",
+    "mae",
+    "mape",
+    "me",
+    "median_absolute_error",
+    "mse",
     "pearson_r",
-    "pearson_r_p_value",
     "pearson_r_eff_p_value",
-    "me",
+    "pearson_r_p_value",
+    "r2",
     "rmse",
-    "mse",
-    "mae",
-    "median_absolute_error",
     "smape",
-    "mape",
     "spearman_r",
-    "spearman_r_p_value",
     "spearman_r_eff_p_value",
-    "effective_sample_size",
-    "r2",
+    "spearman_r_p_value",
 ]
 
 
@@ -71,6 +73,75 @@ def _determine_input_core_dims(dim, weights):
     return input_core_dims
 
 
+def linslope(a, b, dim=None, weights=None, skipna=False, keep_attrs=False):
+    """Slope of linear fit.
+
+    .. math::
+        s_{ab} = \\frac{ \\sum_{i=i}^{n} (a_{i} - \\bar{a}) (b_{i} - \\bar{b}) }
+                 { \\sum_{i=1}^{n} (a_{i} - \\bar{a})^{2} }
+
+    Parameters
+    ----------
+    a : xarray.Dataset or xarray.DataArray
+        Labeled array(s) over which to apply the function.
+    b : xarray.Dataset or xarray.DataArray
+        Labeled array(s) over which to apply the function.
+    dim : str, list
+        The dimension(s) to apply the correlation along. Note that this dimension will
+        be reduced as a result. Defaults to None reducing all dimensions.
+    weights : xarray.Dataset or xarray.DataArray or None
+        Weights matching dimensions of ``dim`` to apply during the function.
+    skipna : bool
+        If True, skip NaNs when computing function.
+    keep_attrs : bool
+        If True, the attributes (attrs) will be copied
+        from the first input to the new one.
+        If False (default), the new object will
+        be returned without attributes.
+
+    Returns
+    -------
+    xarray.DataArray or xarray.Dataset
+        Slope of linear fit.
+
+    See Also
+    --------
+    scipy.stats.linregress
+
+    Examples
+    --------
+    >>> a = xr.DataArray(np.random.rand(5, 3, 3),
+    ...                  dims=['time', 'x', 'y'])
+    >>> b = xr.DataArray(np.random.rand(5, 3, 3),
+    ...                  dims=['time', 'x', 'y'])
+    >>> xs.linslope(a, b, dim='time')
+    <xarray.DataArray (x: 3, y: 3)>
+    array([[-0.30948771, -0.21562529, -0.63141304],
+           [ 0.31446077,  2.23858011,  0.44743617],
+           [-0.22243944,  0.47034784,  1.08512859]])
+    Dimensions without coordinates: x, y
+    """
+    _fail_if_dim_empty(dim)
+    dim, _ = _preprocess_dims(dim, a)
+    a, b = xr.broadcast(a, b, exclude=dim)
+    a, b, new_dim, weights = _stack_input_if_needed(a, b, dim, weights)
+    weights = _preprocess_weights(a, dim, new_dim, weights)
+
+    input_core_dims = _determine_input_core_dims(new_dim, weights)
+
+    return xr.apply_ufunc(
+        _linslope,
+        a,
+        b,
+        weights,
+        input_core_dims=input_core_dims,
+        kwargs={"axis": -1, "skipna": skipna},
+        dask="parallelized",
+        output_dtypes=[float],
+        keep_attrs=keep_attrs,
+    )
+
+
 def pearson_r(a, b, dim=None, weights=None, skipna=False, keep_attrs=False):
     """Pearson's correlation coefficient.
 
diff --git a/xskillscore/core/np_deterministic.py b/xskillscore/core/np_deterministic.py
index 180cbd37..aada365f 100644
--- a/xskillscore/core/np_deterministic.py
+++ b/xskillscore/core/np_deterministic.py
@@ -6,21 +6,22 @@
 from .utils import suppress_warnings
 
 __all__ = [
+    "_effective_sample_size",
+    "_linslope",
+    "_mae",
+    "_mape",
+    "_me",
+    "_median_absolute_error",
+    "_mse",
     "_pearson_r",
-    "_pearson_r_p_value",
     "_pearson_r_eff_p_value",
-    "_me",
+    "_pearson_r_p_value",
+    "_r2",
     "_rmse",
-    "_mse",
-    "_mae",
-    "_median_absolute_error",
     "_smape",
-    "_mape",
     "_spearman_r",
-    "_spearman_r_p_value",
     "_spearman_r_eff_p_value",
-    "_effective_sample_size",
-    "_r2",
+    "_spearman_r_p_value",
 ]
 
 
@@ -153,6 +154,55 @@ def _effective_sample_size(a, b, axis, skipna):
     return n_eff
 
 
+def _linslope(a, b, weights, axis, skipna):
+    """ndarray implementation of scipy.stats.linregress[slope].
+
+    Parameters
+    ----------
+    a : ndarray
+        Input array.
+    b : ndarray
+        Input array.
+    axis : int
+        The axis to apply the linear slope along.
+    weights : ndarray
+        Input array of weights for a and b.
+    skipna : bool
+        If True, skip NaNs when computing function.
+
+    Returns
+    -------
+    res : ndarray
+        slope of linear fit.
+
+    See Also
+    --------
+    scipy.stats.linregress
+    """
+    sumfunc, meanfunc = _get_numpy_funcs(skipna)
+    if skipna:
+        a, b, weights = _match_nans(a, b, weights)
+    weights = _check_weights(weights)
+    a = np.rollaxis(a, axis)
+    b = np.rollaxis(b, axis)
+    if weights is not None:
+        weights = np.rollaxis(weights, axis)
+
+    am, bm = __compute_anomalies(a, b, weights=weights, axis=0, skipna=skipna)
+
+    if weights is not None:
+        s_num = sumfunc(weights * am * bm, axis=0)
+        s_den = sumfunc(weights * am * am, axis=0)
+    else:
+        s_num = sumfunc(am * bm, axis=0)
+        s_den = sumfunc(am * am, axis=0)
+
+    with suppress_warnings("invalid value encountered in true_divide"):
+        with suppress_warnings("invalid value encountered in double_scalars"):
+            res = s_num / s_den
+    return res
+
+
 def _r2(a, b, weights, axis, skipna):
     """ndarray implementation of sklearn.metrics.r2_score.
 
diff --git a/xskillscore/tests/test_accessor_deterministic.py b/xskillscore/tests/test_accessor_deterministic.py
index 1dd0b7bd..d36271b8 100644
--- a/xskillscore/tests/test_accessor_deterministic.py
+++ b/xskillscore/tests/test_accessor_deterministic.py
@@ -4,6 +4,7 @@
 
 from xskillscore.core.deterministic import (
     effective_sample_size,
+    linslope,
     mae,
     mape,
     me,
@@ -21,6 +22,7 @@
 )
 
 correlation_metrics = [
+    linslope,
     pearson_r,
     r2,
     pearson_r_p_value,
diff --git a/xskillscore/tests/test_deterministic.py b/xskillscore/tests/test_deterministic.py
index 0becc45c..bf97a1ae 100644
--- a/xskillscore/tests/test_deterministic.py
+++ b/xskillscore/tests/test_deterministic.py
@@ -7,6 +7,7 @@
     _preprocess_dims,
     _preprocess_weights,
     effective_sample_size,
+    linslope,
     mae,
     mape,
     me,
@@ -24,6 +25,7 @@
 )
 from xskillscore.core.np_deterministic import (
     _effective_sample_size,
+    _linslope,
     _mae,
     _mape,
     _me,
@@ -41,6 +43,7 @@
 )
 
 correlation_metrics = [
+    (linslope, _linslope),
     (pearson_r, _pearson_r),
     (r2, _r2),
     (pearson_r_p_value, _pearson_r_p_value),
@@ -51,6 +54,7 @@
     (effective_sample_size, _effective_sample_size),
 ]
 correlation_metrics_names = [
+    "linslope",
     "pearson_r",
     "r2",
     "pearson_r_p_value",
diff --git a/xskillscore/tests/test_gridded_metrics.py b/xskillscore/tests/test_gridded_metrics.py
index 180abe14..cf9009fb 100644
--- a/xskillscore/tests/test_gridded_metrics.py
+++ b/xskillscore/tests/test_gridded_metrics.py
@@ -2,6 +2,7 @@
 import pytest
 
 from xskillscore.core.deterministic import (
+    linslope,
     mae,
     mape,
     me,
@@ -17,6 +18,7 @@
 )
 
 METRICS = [
+    linslope,
     mae,
     mse,
     median_absolute_error,
diff --git a/xskillscore/tests/test_mask_skipna.py b/xskillscore/tests/test_mask_skipna.py
index 879e3b90..7f038e4e 100644
--- a/xskillscore/tests/test_mask_skipna.py
+++ b/xskillscore/tests/test_mask_skipna.py
@@ -2,6 +2,7 @@
 import pytest
 
 from xskillscore.core.deterministic import (
+    linslope,
     mae,
     mape,
     median_absolute_error,
diff --git a/xskillscore/tests/test_metric_results_accurate.py b/xskillscore/tests/test_metric_results_accurate.py
index 88995408..8b402099 100644
--- a/xskillscore/tests/test_metric_results_accurate.py
+++ b/xskillscore/tests/test_metric_results_accurate.py
@@ -1,7 +1,7 @@
 import numpy as np
 import pytest
 import sklearn.metrics
-from scipy.stats import pearsonr, spearmanr
+from scipy.stats import linregress, pearsonr, spearmanr
 from sklearn.metrics import (
     mean_absolute_error,
     mean_absolute_percentage_error,
@@ -11,6 +11,7 @@
 
 import xskillscore as xs
 from xskillscore.core.deterministic import (
+    linslope,
     mae,
     mape,
     me,
@@ -36,6 +37,7 @@
 ]
 
 xs_scipy_metrics = [
+    (linslope, linregress, 0),
     (pearson_r, pearsonr, 0),
     (spearman_r, spearmanr, 0),
     (pearson_r_p_value, pearsonr, 1),
diff --git a/xskillscore/tests/test_skipna_functionality.py b/xskillscore/tests/test_skipna_functionality.py
index 50adb726..7024cf9d 100644
--- a/xskillscore/tests/test_skipna_functionality.py
+++ b/xskillscore/tests/test_skipna_functionality.py
@@ -4,6 +4,7 @@
 from xarray.tests import CountingScheduler, assert_allclose, raise_if_dask_computes
 
 from xskillscore.core.deterministic import (
+    linslope,
     mae,
     mape,
     me,
@@ -19,6 +20,7 @@
 )
 
 WEIGHTED_METRICS = [
+    linslope,
     pearson_r,
     pearson_r_p_value,
     spearman_r,
diff --git a/xskillscore/tests/test_weighted_metric_results_accurate.py b/xskillscore/tests/test_weighted_metric_results_accurate.py
index 5ea8e637..6de2338e 100644
--- a/xskillscore/tests/test_weighted_metric_results_accurate.py
+++ b/xskillscore/tests/test_weighted_metric_results_accurate.py
@@ -1,7 +1,5 @@
 import numpy as np
 import pytest
-import sklearn.metrics
-from scipy.stats import pearsonr
 from sklearn.metrics import (
     mean_absolute_error,
     mean_absolute_percentage_error,
@@ -11,6 +9,7 @@
 
 import xskillscore as xs
 from xskillscore.core.deterministic import (
+    linslope,
     mae,
     mape,
     me,
@@ -44,7 +43,15 @@ def weighted_pearsonr(x, y, w):
     return r_num / r_den
 
 
-xs_scipy_metrics = [(pearson_r, weighted_pearsonr)]
+def weighted_linslope(x, y, w):
+    xm = x - (np.sum(x * w) / np.sum(w))
+    ym = y - (np.sum(y * w) / np.sum(w))
+    s_num = np.sum(w * xm * ym)
+    s_den = np.sum(w * xm * xm)
+    return s_num / s_den
+
+
+xs_scipy_metrics = [(pearson_r, weighted_pearsonr), (linslope, weighted_linslope)]
 
 
 xs_np_metrics = [