include xr-ed versions of properscorings crps_gaussian, crps_ensemble

README.rst

@@ -4,7 +4,7 @@ xskillscore: Metrics for verifying forecasts
 .. image:: https://travis-ci.org/raybellwaves/xskillscore.svg?branch=master
    :target: https://travis-ci.org/raybellwaves/xskillscore
 
-**xskillscore** is an open source project and Python package that provides verification metrics of deterministic (and probabilistic) forecasts with xarray.
+**xskillscore** is an open source project and Python package that provides verification metrics of deterministic (and probabilistic from `properscoring`) forecasts with `xarray`.
 
 Installing
 ----------
@@ -32,9 +32,10 @@ Examples
                       dims=['time', 'lat', 'lon'])
    fct = xr.DataArray(np.random.rand(3, 4, 5),
                       coords=[pd.date_range('1/1/2000', '1/3/2000', freq='D'),
-                              np.arange(4), np.arange(5)],         
+                              np.arange(4), np.arange(5)],
                       dims=['time', 'lat', 'lon'])
-
+
+   # deterministic
    r = xs.pearson_r(obs, fct, 'time')
    #>>> r
    #<xarray.DataArray (lat: 4, lon: 5)>
@@ -45,11 +46,19 @@ Examples
    #Coordinates:
    #  * lat      (lat) int64 0 1 2 3
    #  * lon      (lon) int64 0 1 2 3 4
-  
+
    r_p_value = xs.pearson_r_p_value(obs, fct, 'time')
 
    rmse = xs.rmse(obs, fct, 'time')
 
    mse = xs.mse(obs, fct, 'time')
 
-   mae = xs.mae(obs, fct, 'time') 
+   mae = xs.mae(obs, fct, 'time')
+
+
+   # probabilistic
+   crps_ensemble = xs.crps_ensemble(obs, fct)
+
+   crps_gaussian = xs.crps_gaussian(obs, fct.mean('time'), fct.std('time'))
+
+   threshold_brier_score = xs.threshold_brier_score(obs, fct, 0.7)

ci/requirements-py36.yml

@@ -8,3 +8,5 @@ dependencies:
   - dask
   - scipy
   - pytest
+  - pip:
+    - properscoring

setup.py

@@ -1,6 +1,5 @@
 from setuptools import find_packages, setup
 
-
 DISTNAME = 'xskillscore'
 VERSION = '0.0.2'
 LICENSE = 'Apache'
@@ -9,7 +8,7 @@
 DESCRIPTION = "xskillscore"
 LONG_DESCRIPTION = """Metrics for verifying forecasts"""
 URL = 'https://github.com/raybellwaves/xskillscore'
-INSTALL_REQUIRES = ['scikit-learn', 'xarray', 'dask', 'scipy']
+INSTALL_REQUIRES = ['scikit-learn', 'xarray', 'dask', 'scipy', 'properscoring']
 TESTS_REQUIRE = ['pytest']
 PYTHON_REQUIRE = '>=3.6'
 

xskillscore/__init__.py

@@ -1,2 +1,6 @@
 # flake8: noqa
-from .core.deterministic import pearson_r, pearson_r_p_value, rmse, mse, mae
+from .core.deterministic import mae, mse, pearson_r, pearson_r_p_value, rmse
+from .core.probabilistic import xr_crps_ensemble as crps_ensemble
+from .core.probabilistic import xr_crps_gaussian as crps_gaussian
+from .core.probabilistic import \
+    xr_threshold_brier_score as threshold_brier_score

xskillscore/core/probabilistic.py

@@ -0,0 +1,129 @@
+import xarray as xr
+from properscoring import crps_ensemble, crps_gaussian, threshold_brier_score
+
+
+def xr_crps_gaussian(observations, mu, sig):
+    """
+    xarray version of properscoring.crps_gaussian.
+
+    Parameters
+    ----------
+    observations : Dataset, DataArray, GroupBy, Variable, numpy/dask arrays or
+     scalars, Mix of labeled and/or unlabeled observations arrays.
+    mu : Dataset, DataArray, GroupBy, Variable, numpy/dask arrays or
+     scalars, Mix of labeled and/or unlabeled forecasts mean arrays.
+    sig : Dataset, DataArray, GroupBy, Variable, numpy/dask arrays or
+     scalars, Mix of labeled and/or unlabeled forecasts mean arrays.
+
+    Returns
+    -------
+    Single value or tuple of Dataset, DataArray, Variable, dask.array.Array or
+     numpy.ndarray, the first type on that list to appear on an input.
+
+    See Also
+    --------
+    properscoring.crps_gaussian
+    xarray.apply_ufunc
+    """
+    # check if same dimensions
+    if mu.dims != observations.dims:
+        observations, mu = xr.broadcast(observations, mu)
+    if sig.dims != observations.dims:
+        observations, sig = xr.broadcast(observations, sig)
+    return xr.apply_ufunc(crps_gaussian,
+                          observations,
+                          mu,
+                          sig,
+                          input_core_dims=[[], [], []],
+                          dask='parallelized',
+                          output_dtypes=[float])
+
+
+def xr_crps_ensemble(observations, forecasts):
+    """
+    xarray version of properscoring.crps_ensemble.
+
+    Parameters
+    ----------
+    observations : Dataset, DataArray, GroupBy, Variable, numpy/dask arrays or
+     scalars, Mix of labeled and/or unlabeled observations arrays.
+    forecasts : Dataset, DataArray, GroupBy, Variable, numpy/dask arrays or
+     scalars, Mix of labeled and/or unlabeled forecasts arrays.
+
+    Returns
+    -------
+    Single value or tuple of Dataset, DataArray, Variable, dask.array.Array or
+    numpy.ndarray, the first type on that list to appear on an input.
+
+    See Also
+    --------
+    properscoring.crps_ensemble
+    xarray.apply_ufunc
+    """
+    if forecasts.dims != observations.dims:
+        observations, forecasts = xr.broadcast(observations, forecasts)
+    return xr.apply_ufunc(crps_ensemble,
+                          observations,
+                          forecasts,
+                          input_core_dims=[[], []],
+                          dask='parallelized',
+                          output_dtypes=[float])
+
+
+def xr_threshold_brier_score(observations,
+                             forecasts,
+                             threshold,
+                             issorted=False,
+                             axis=-1):
+    """
+    xarray version of properscoring.threshold_brier_score: Calculate the Brier
+     scores of an ensemble for exceeding given thresholds.
+
+    Parameters
+    ----------
+    observations : Dataset, DataArray, GroupBy, Variable, numpy/dask arrays or
+     scalars, Mix of labeled and/or unlabeled observations arrays.
+    forecasts : Dataset, DataArray, GroupBy, Variable, numpy/dask arrays or
+     scalars, Mix of labeled and/or unlabeled forecasts arrays.
+    threshold : scalar (not yet implemented: or 1d scalar threshold value(s) at
+     which to calculate) exceedence Brier scores.
+    issorted : bool, optional
+        Optimization flag to indicate that the elements of `ensemble` are
+        already sorted along `axis`.
+    axis : int, optional
+        Axis in forecasts which corresponds to different ensemble members,
+        along which to calculate the threshold decomposition.
+
+
+    Returns
+    -------
+    Single value or tuple of Dataset, DataArray, Variable, dask.array.Array or
+    numpy.ndarray, the first type on that list to appear on an input. (If
+    ``threshold`` is a scalar, the result will have the same shape as
+    observations. Otherwise, it will have an additional final dimension
+    corresponding to the threshold levels.)
+
+    References
+    ----------
+    Gneiting, T. and Ranjan, R. Comparing density forecasts using threshold-
+       and quantile-weighted scoring rules. J. Bus. Econ. Stat. 29, 411-422
+       (2011). http://www.stat.washington.edu/research/reports/2008/tr533.pdf
+
+    See Also
+    --------
+    properscoring.threshold_brier_score
+    xarray.apply_ufunc
+    """
+    if forecasts.dims != observations.dims:
+        observations, forecasts = xr.broadcast(observations, forecasts)
+    return xr.apply_ufunc(threshold_brier_score,
+                          observations,
+                          forecasts,
+                          threshold,
+                          input_core_dims=[[], [], []],
+                          kwargs={
+                              'axis': axis,
+                              'issorted': issorted
+                          },
+                          dask='parallelized',
+                          output_dtypes=[float])

xskillscore/tests/test_probabilistic.py

@@ -0,0 +1,71 @@
+import numpy as np
+import pandas as pd
+import pytest
+import xarray as xr
+from properscoring import (crps_ensemble, crps_gaussian,
+                           threshold_brier_score)
+from xarray.tests import assert_identical
+
+from xskillscore.core.probabilistic import (xr_crps_ensemble, xr_crps_gaussian,
+                                            xr_threshold_brier_score)
+
+
+@pytest.fixture
+def a_dask():
+    dates = pd.date_range('1/1/2000', '1/3/2000', freq='D')
+    lats = np.arange(4)
+    lons = np.arange(5)
+    data = np.random.rand(len(dates), len(lats), len(lons))
+    return xr.DataArray(data,
+                        coords=[dates, lats, lons],
+                        dims=['time', 'lat', 'lon']).chunk()
+
+
+@pytest.fixture
+def b_dask():
+    dates = pd.date_range('1/1/2000', '1/3/2000', freq='D')
+    lats = np.arange(4)
+    lons = np.arange(5)
+    data = np.random.rand(len(dates), len(lats), len(lons))
+    return xr.DataArray(data,
+                        coords=[dates, lats, lons],
+                        dims=['time', 'lat', 'lon']).chunk()
+
+
+def test_xr_crps_ensemble_dask(a_dask, b_dask):
+    actual = xr_crps_ensemble(a_dask, b_dask)
+    expected = crps_ensemble(a_dask, b_dask)
+    expected = xr.DataArray(expected, coords=a_dask.coords)
+    # test for numerical identity of xr_crps and crps
+    assert_identical(actual, expected)
+    # test that xr_crps_ensemble returns chunks
+    assert actual.chunks is not None
+    # show that crps_ensemble returns no chunks
+    assert expected.chunks is None
+
+
+def test_xr_crps_gaussian_dask(a_dask, b_dask):
+    mu = b_dask.mean('time')
+    sig = b_dask.std('time')
+    actual = xr_crps_gaussian(a_dask, mu, sig)
+    expected = crps_gaussian(a_dask, mu, sig)
+    expected = xr.DataArray(expected, coords=a_dask.coords)
+    # test for numerical identity of xr_crps and crps
+    assert_identical(actual, expected)
+    # test that xr_crps_ensemble returns chunks
+    assert actual.chunks is not None
+    # show that crps_ensemble returns no chunks
+    assert expected.chunks is None
+
+
+def test_xr_threshold_brier_score_dask(a_dask, b_dask):
+    threshold = .5
+    actual = xr_threshold_brier_score(a_dask, b_dask, threshold)
+    expected = threshold_brier_score(a_dask, b_dask, threshold)
+    expected = xr.DataArray(expected, coords=a_dask.coords)
+    # test for numerical identity of xr_threshold and threshold
+    assert_identical(actual, expected)
+    # test that xr_crps_ensemble returns chunks
+    assert actual.chunks is not None
+    # show that crps_ensemble returns no chunks
+    assert expected.chunks is None