Repository: Ouranosinc/xclim
Pull request #29 (2018-09-25 14:06:45)

initial try to implement the prcp_tot indice

coveralls: ## Pull Request Test Coverage Report for [Build 88](https://coveralls.io/builds/19179396)

* **14** of **25**   **(56.0%)**  changed or added relevant lines in **3** files are covered.
* No unchanged relevant lines lost coverage.
* Overall coverage increased (+**0.08%**) to **51.957%**

---

|  Changes Missing Coverage | Covered Lines | Changed/Added Lines | % |
| :-----|--------------|--------|---: |
| [xclim/checks.py](https://coveralls.io/builds/19179396/source?filename=xclim%2Fchecks.py#L110) | 2 | 7 | 28.57%
| [xclim/indices.py](https://coveralls.io/builds/19179396/source?filename=xclim%2Findices.py#L69) | 11 | 17 | 64.71%
<!-- | **Total:** | **14** | **25** | **56.0%** | -->


|  Totals | [![Coverage Status](https://coveralls.io/builds/19179396/badge)](https://coveralls.io/builds/19179396) |
| :-- | --: |
| Change from base [Build 85](https://coveralls.io/builds/19127714): |  0.08% |
| Covered Lines: | 146 |
| Relevant Lines: | 281 |

---
##### 💛  - [Coveralls](https://coveralls.io)

codecov[bot]: # [Codecov](https://codecov.io/gh/Ouranosinc/xclim/pull/29?src=pr&el=h1) Report
> Merging [#29](https://codecov.io/gh/Ouranosinc/xclim/pull/29?src=pr&el=desc) into [master](https://codecov.io/gh/Ouranosinc/xclim/commit/2ecc7a852b1abc233248c88ea895694fe8d97bfb?src=pr&el=desc) will **increase** coverage by `0.07%`.
> The diff coverage is `56%`.

[![Impacted file tree graph](https://codecov.io/gh/Ouranosinc/xclim/pull/29/graphs/tree.svg?width=650&token=1clOxQaUPR&height=150&src=pr)](https://codecov.io/gh/Ouranosinc/xclim/pull/29?src=pr&el=tree)

```diff
@@            Coverage Diff             @@
##           master      #29      +/-   ##
==========================================
+ Coverage   51.87%   51.95%   +0.07%     
==========================================
  Files           5        5              
  Lines         266      281      +15     
==========================================
+ Hits          138      146       +8     
- Misses        128      135       +7
```


| [Impacted Files](https://codecov.io/gh/Ouranosinc/xclim/pull/29?src=pr&el=tree) | Coverage Δ | |
|---|---|---|
| [xclim/\_\_init\_\_.py](https://codecov.io/gh/Ouranosinc/xclim/pull/29/diff?src=pr&el=tree#diff-eGNsaW0vX19pbml0X18ucHk=) | `76.92% <100%> (ø)` | :arrow_up: |
| [xclim/checks.py](https://codecov.io/gh/Ouranosinc/xclim/pull/29/diff?src=pr&el=tree#diff-eGNsaW0vY2hlY2tzLnB5) | `64.15% <28.57%> (-6.07%)` | :arrow_down: |
| [xclim/indices.py](https://codecov.io/gh/Ouranosinc/xclim/pull/29/diff?src=pr&el=tree#diff-eGNsaW0vaW5kaWNlcy5weQ==) | `55.69% <64.7%> (+1.33%)` | :arrow_up: |

------

[Continue to review full report at Codecov](https://codecov.io/gh/Ouranosinc/xclim/pull/29?src=pr&el=continue).
> **Legend** - [Click here to learn more](https://docs.codecov.io/docs/codecov-delta)
> `Δ = absolute <relative> (impact)`, `ø = not affected`, `? = missing data`
> Powered by [Codecov](https://codecov.io/gh/Ouranosinc/xclim/pull/29?src=pr&el=footer). Last update [2ecc7a8...30f97ab](https://codecov.io/gh/Ouranosinc/xclim/pull/29?src=pr&el=lastupdated). Read the [comment docs](https://docs.codecov.io/docs/pull-request-comments).

Balinus: I suggest `prcp_tot` to keep inline with xclim usage of naming functions. e.g. see `summer_days`

Patch:
diff --git a/.travis.yml b/.travis.yml
index 09f889cc..9a6116f1 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -29,10 +29,11 @@ matrix:
       python: 2.7
 
 # Command to install dependencies, e.g. pip install -r requirements.txt --use-mirrors
-install: pip install -U tox-travis
+install: 
+  - pip install -U tox-travis 
 
 # Command to run tests, e.g. python setup.py test
-script: tox
+script: travis_wait tox -e $TOXENV
 
 branches:
   only:
diff --git a/setup.cfg b/setup.cfg
index f6788a0c..0768ba68 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 0.4-alpha
+current_version = 0.5-alpha
 commit = True
 tag = True
 parse = (?P<major>\d+)\.(?P<minor>\d+)(\-(?P<release>[a-z]+))?
diff --git a/setup.py b/setup.py
index 8bca66ec..b163960b 100644
--- a/setup.py
+++ b/setup.py
@@ -54,7 +54,7 @@ setup(
     tests_require=test_requirements,
     extras_require={'docs': docs_requirements},
     url='https://github.com/Ouranosinc/xclim',
-    version='0.4-alpha',
+    version='0.5-alpha',
     zip_safe=False,
 )
 
diff --git a/xclim/__init__.py b/xclim/__init__.py
index 61e06a36..2e67b9ee 100644
--- a/xclim/__init__.py
+++ b/xclim/__init__.py
@@ -4,7 +4,7 @@
 
 __author__ = """Travis Logan"""
 __email__ = 'logan.travis@ouranos.ca'
-__version__ = '0.4-alpha'
+__version__ = '0.5-alpha'
 
 from . import indices
 # from .stats import fit, test
diff --git a/xclim/checks.py b/xclim/checks.py
index a8c964ed..fe29bbc3 100644
--- a/xclim/checks.py
+++ b/xclim/checks.py
@@ -1,5 +1,5 @@
-from warnings import warn
 from functools import wraps
+from warnings import warn
 
 """
 Dev notes
@@ -101,3 +101,16 @@ def valid_missing_data_threshold(comp, threshold=0):
     """Check that the relative number of missing data points does not exceed a threshold."""
     # TODO
     pass
+
+
+def check_is_dataarray(comp):
+    """Decorator to check that a computation has an instance of xarray.DataArray
+     as first argument."""
+
+    @wraps(comp)
+    def func(data_array, *args, **kwds):
+        import xarray as xr
+        assert (isinstance(data_array, xr.DataArray))
+        return comp(data_array, *args, **kwds)
+
+    return func
diff --git a/xclim/indices.py b/xclim/indices.py
index d2e6d8cb..d1d27e03 100644
--- a/xclim/indices.py
+++ b/xclim/indices.py
@@ -2,19 +2,21 @@
 
 """Main module.
 """
+from functools import wraps
+
 import numpy as np
 import xarray as xr
 
+from . import run_length as rl
 from .checks import valid_daily_mean_temperature, valid_daily_max_min_temperature, valid_daily_min_temperature, \
     valid_daily_max_temperature, valid_daily_mean_discharge
-from . import run_length as rl
-from functools import wraps
 
 # Frequencies : YS: year start, QS-DEC: seasons starting in december, MS: month start
 # See http://pandas.pydata.org/pandas-docs/stable/timeseries.html#offset-aliases
 K2C = 273.15
 ftomm = np.nan
 
+
 # TODO: Move utility functions to another file.
 
 
@@ -32,6 +34,7 @@ def with_attrs(**func_attrs):
     -----
     Assumes the output has an attrs dictionary attribute (e.g. xarray.DataArray).
     """
+
     def attr_decorator(fn):
         # Use the docstring as the description attribute.
         func_attrs['description'] = first_paragraph(fn.__doc__)
@@ -50,6 +53,7 @@ def with_attrs(**func_attrs):
 
     return attr_decorator
 
+
 # -------------------------------------------------- #
 # ATTENTION: ASSUME ALL INDICES WRONG UNTIL TESTED ! #
 # -------------------------------------------------- #
@@ -62,7 +66,7 @@ def base_flow_index(q, freq='YS'):
     mq = q.resample(time=freq)
 
     m7m = m7.min(dim='time')
-    return m7m/mq.mean(dim='time')
+    return m7m / mq.mean(dim='time')
 
 
 @valid_daily_min_temperature
@@ -89,6 +93,7 @@ def cold_spell_duration_index(tasmin, tn10, freq='YS'):
     --------
     percentile_doy
     """
+
     def func(x):
         xr.apply_ufunc(rl.windowed_run_count,
                        x,
@@ -224,9 +229,9 @@ def CFD2(tasmin, freq='AS-JUL'):
 def cooling_degree_days(tas, thresh=18, freq='YS'):
     """Cooling degree days above threshold."""
 
-    return tas.pipe(lambda x: x - thresh - K2C)\
-        .clip(min=0)\
-        .resample(time=freq)\
+    return tas.pipe(lambda x: x - thresh - K2C) \
+        .clip(min=0) \
+        .resample(time=freq) \
         .sum(dim='time')
 
 
@@ -268,14 +273,14 @@ def frost_days(tasmin, freq='YS'):
 
 
 @valid_daily_mean_temperature
-def growing_degree_days(tas, thresh=4,  freq='YS'):
+def growing_degree_days(tas, thresh=4, freq='YS'):
     """Growing degree days over 4℃.
 
     The sum of degree-days over 4℃.
     """
-    return tas.pipe(lambda x: x - thresh - K2C)\
-        .clip(min=0)\
-        .resample(time=freq)\
+    return tas.pipe(lambda x: x - thresh - K2C) \
+        .clip(min=0) \
+        .resample(time=freq) \
         .sum(dim='time')
 
 
@@ -353,9 +358,9 @@ def heating_degree_days(tas, freq='YS', thresh=17):
 
     Sum of positive values of threshold - tas.
     """
-    return tas.pipe(lambda x: thresh - x)\
-        .clip(0)\
-        .resample(time=freq)\
+    return tas.pipe(lambda x: thresh - x) \
+        .clip(0) \
+        .resample(time=freq) \
         .sum(dim='time')
 
 
@@ -449,11 +454,70 @@ def tn_min(tasmin, freq='YS'):
     return tasmin.resample(time=freq).min(dim='time')
 
 
+# @check_is_dataarray
+def prcp_tot(pr, freq='YS', units='kg m-2 s-1'):
+    r"""Accumulated total (liquid + solid) precipitation
+
+    Resample the original daily mean precipitation flux and cumulate over each period
+
+    Parameters
+    ----------
+    pr : xarray.DataArray
+      Mean daily precipitation flux [Kg m-2 s-1] or [mm].
+    freq : str, optional
+      Resampling frequency as defined in
+      http://pandas.pydata.org/pandas-docs/stable/timeseries.html#resampling.
+    units: str, optional
+      units of the precipitation data. Must be within ['kg m-2 s-2', 'mm']
+
+    Returns
+    -------
+    xarray.DataArray
+      The total daily precipitation at the given time frequency in [mm].
+
+
+    Notes
+    -----
+    Let :math:`T_i` be the mean daily temperature of day `i`, then for a period `p` starting at
+    day `a` and finishing on day `b`
+
+    .. math::
+
+       TG_p = \frac{\sum_{i=a}^{b} T_i}{b - a + 1}
+
+
+    Examples
+    --------
+    The following would compute for each grid cell of file `pr_day.nc` the total
+    precipitation at the seasonal frequency, ie DJF, MAM, JJA, SON, DJF, etc.
+
+    >>> pr_day = xr.open_dataset('pr_day.nc').pr
+    >>> prcp_tot_seasonal = prcp_tot(pr_day, freq="QS-DEC")
+
+    """
+    # TODO deal with the time_boundaries
+
+    # resample the precipitation to the wanted frequency
+    arr = pr.resample(time=freq)
+    # cumulate the values over the season
+    output = arr.sum(dim='time')
+    # unit conversion as needed
+    if units == 'kg m-2 s-1':
+        # convert from km m-2 s-1 to mm day-1
+        output *= 86400  # number of sec in 24h
+    elif units == 'mm':
+        # nothing to do
+        pass
+    else:
+        raise RuntimeError('non conform units')
+    return output
+
+
 @valid_daily_min_temperature
 def tropical_nights(tasmin, thresh=20, freq='YS'):
     """Number of days with minimum daily temperature above threshold."""
-    return tasmin.pipe(lambda x: (tasmin > thresh + K2C) * 1)\
-        .resample(time=freq)\
+    return tasmin.pipe(lambda x: (tasmin > thresh + K2C) * 1) \
+        .resample(time=freq) \
         .sum(dim='time')
 
 


Test patch:
diff --git a/tests/test_indices.py b/tests/test_indices.py
index dc6917dd..45310ff3 100644
--- a/tests/test_indices.py
+++ b/tests/test_indices.py
@@ -17,6 +17,8 @@ the results in a reference netcdf dataset. We could then compare the hailstorm o
 of defense.
 """
 
+# import cftime
+import calendar
 import os
 
 import numpy as np
@@ -84,6 +86,28 @@ class Test_cooling_degree_days():
         assert len(cdd.description) > 0
 
 
+class Test_prcptotal():
+    # build test data for different calendar
+    time_std = pd.date_range('2000-01-01', '2010-12-31', freq='D')
+    da_std = xr.DataArray(time_std.year, coords=[time_std], dims='time')
+
+    # calendar 365_day and 360_day not tested for now since xarray.resample
+    # does not support other calendars than standard
+    #
+    # units = 'days since 2000-01-01 00:00'
+    # time_365 = cftime.num2date(np.arange(0, 10 * 365), units, '365_day')
+    # time_360 = cftime.num2date(np.arange(0, 10 * 360), units, '360_day')
+    # da_365 = xr.DataArray(np.arange(time_365.size), coords=[time_365], dims='time')
+    # da_360 = xr.DataArray(np.arange(time_360.size), coords=[time_360], dims='time')
+
+    def test_yearly(self):
+        da_std = self.da_std
+        out_std = xci.prcp_tot(da_std, units='mm')
+        # l_years = np.unique(da_std.time.dt.year) TODO: Unused local variables are a PEP8 violation
+        target = [(365 + calendar.isleap(y)) * y for y in np.unique(da_std.time.dt.year)]
+        assert (np.allclose(target, out_std.values))
+
+
 # I'd like to parametrize some of these tests so we don't have to write individual tests for each indicator.
 @pytest.mark.skip('')
 class TestTG():