Repository: Ouranosinc/xclim

Issue:
Indicator for average temperature based on tasmin and tasmax
### Description
<!--Describe what you were trying to get done.
Tell us what happened, what went wrong, and what you expected to happen.-->
Create a new indicator computing the daily average temperature from the daily min and max temperature. 


Discussion:
coveralls: ## Pull Request Test Coverage Report for [Build 1278](https://coveralls.io/builds/27259303)

* **3** of **3**   **(100.0%)**  changed or added relevant lines in **1** file are covered.
* No unchanged relevant lines lost coverage.
* Overall coverage increased (+**0.02%**) to **90.786%**

---



|  Totals | [![Coverage Status](https://coveralls.io/builds/27259303/badge)](https://coveralls.io/builds/27259303) |
| :-- | --: |
| Change from base [Build 1260](https://coveralls.io/builds/27178478): |  0.02% |
| Covered Lines: | 3094 |
| Relevant Lines: | 3408 |

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

aulemahal: Added the atmos indicator. However, I add to hack the `Indicator.missing()` method. Right now it only checks for missing days in the data. But that doesn't make sense for not-resampling indicators, like this one. So I return "no missing values" if "freq" is not an argument. This could be modified if we add non-resampling missing value checks later on.
huard: Another option would be to subclass the missing method with a function that returns False. 
That being said, I think we should add a test to check that when a tasmin is missing but not a tasmax, the result is a NaN. 
aulemahal: 1) You mean like:
```python
class TasminTasmaxNoResample(TasminTasmax):
    @staticmethod
    def missing *args, **kwds):
        return False
```
? Seems a very specific solver for a more general problem, no?

Also, xarray already does the check for an addition:
```python
a = xr.DataArray([1, np.nan, 1], dims=('x',), coords={'x': [0, 1, 2]})
b = xr.DataArray([1, 2, 1, 1], dims=('x',), coords={'x': [0, 1, 2, 3]})
a + b
```
gives: `xarray.DataArray (x: 3)> array([ 2., nan,  2.])  Coordinates:   * x        (x) int64 0 1 2`
So, missing values are omitted and NaNs are propagated.
huard: Yes exactly. To me this communicates better what is going on to readers. Happy to be convinced otherwise. @tlvu Any opinion on this ? 

Perfect ! Just want to make sure that as these inputs go through all the Indicator methods, the same is true. They should in theory, but in practice theory is different from practice. 
tlvu: > Another option would be to subclass the missing method with a function that returns False.
> That being said, I think we should add a test to check that when a tasmin is missing but not a tasmax, the result is a NaN.
> Yes exactly. To me this communicates better what is going on to readers. Happy to be convinced otherwise. @tlvu Any opinion on this ?

Having a static method that always return False is weird to me but if that's how the API/library behave/expect then fine by me.

tlogan2000: > To be used as a WPS service, we have to expose an indicator in atmos. I think that was what Travis was looking for. Can be another PR.

Yes. The problem for many of our datasets is that we often estimate 'tas' from 'tasmin' and 'tasmax'.  Using ```xclim``` / ```xarray``` we can do this on the fly and then calculate any index requiring tas as input .  However calling the same calculation via the finch wps server requires the resource netcdf to have the input variable already.

We might want to think about creating/writing to disk the a set of 'tas' .nc files for our own datasets on the thredds server.  However a vision of remotely connecting to other nodes will likely result in a similar issue so this will still be of use I think
aulemahal: I see!
Here I come with another proposition. @huard , it goes to your first comment and as the __add__ of xarray already does this, it's kinda useless for this specific case. But, as long as atmos requires input to be daily, I think giving a default `freq='D'` to the missing check makes sense. And as it flags values missing from the input dataset, still does what we expect it to.
So I am pushing away the issue to the future, when we want to make atmos more general :).

Patch:
diff --git a/HISTORY.rst b/HISTORY.rst
index 4994460e..3ff1c54d 100644
--- a/HISTORY.rst
+++ b/HISTORY.rst
@@ -8,7 +8,7 @@ History
 * Added a `tolerance` argument to `subset_gridpoint` raising an error if distance to closest point is larger than tolerance.
 * Created land module for standardized access to streamflow indices.
 * Enhancement to utils.Indicator to have more dynamic attributes using callables.
-* Added indice `heat_wave_total_length`.
+* Added indices `heat_wave_total_length` and `tas` / `tg`.
 * Fixed a bug with typed call signatures that caused downstream failures on library import.
 
 0.11.x-beta (2019-10-17)
diff --git a/xclim/atmos/_temperature.py b/xclim/atmos/_temperature.py
index 9f645479..642e5094 100644
--- a/xclim/atmos/_temperature.py
+++ b/xclim/atmos/_temperature.py
@@ -13,6 +13,7 @@ __all__ = [
     "heat_wave_frequency",
     "heat_wave_max_length",
     "heat_wave_index",
+    "tg",
     "tg_mean",
     "tg10p",
     "tg90p",
@@ -186,6 +187,16 @@ heat_wave_index = Tasmax(
     compute=indices.heat_wave_index,
 )
 
+tg = TasminTasmax(
+    identifier="tg",
+    units="K",
+    standard_name="air_temperature",
+    long_name="Daily mean temperature",
+    description="Estimated mean temperature from maximum and minimum temperatures",
+    cell_methods="",
+    compute=indices.tas,
+)
+
 tg_mean = Tas(
     identifier="tg_mean",
     units="K",
diff --git a/xclim/indices/_multivariate.py b/xclim/indices/_multivariate.py
index a2671e7f..78d7e145 100644
--- a/xclim/indices/_multivariate.py
+++ b/xclim/indices/_multivariate.py
@@ -33,6 +33,7 @@ __all__ = [
     "liquid_precip_ratio",
     "precip_accumulation",
     "rain_on_frozen_ground_days",
+    "tas",
     "tg90p",
     "tg10p",
     "tn90p",
@@ -801,6 +802,27 @@ def fraction_over_precip_thresh(
     return over / total
 
 
+@declare_units("K", tasmin="[temperature]", tasmax="[temperature]")
+def tas(tasmin: xarray.DataArray, tasmax: xarray.DataArray) -> xarray.DataArray:
+    """Average temperature from minimum and maximum temperatures.
+
+    We assume a symmetrical distribution for the temperature and retrieve the average value as Tg = (Tx + Tn) / 2
+
+    Parameters
+    ----------
+    tasmin : xarray.DataArray
+        Minimum (daily) temperature [℃] or [K]
+    tasmax : xarray.DataArray
+        Maximum (daily) temperature [℃] or [K]
+
+    Returns
+    -------
+    xarray.DataArray
+        Mean (daily) temperature [°C]
+    """
+    return (tasmax + tasmin) / 2
+
+
 @declare_units("days", tas="[temperature]", t90="[temperature]")
 def tg90p(
     tas: xarray.DataArray, t90: xarray.DataArray, freq: str = "YS"
diff --git a/xclim/utils.py b/xclim/utils.py
index 4f9a7f35..17e7f445 100644
--- a/xclim/utils.py
+++ b/xclim/utils.py
@@ -958,7 +958,8 @@ class Indicator:
         """Return whether an output is considered missing or not."""
         from functools import reduce
 
-        freq = kwds.get("freq")
+        freq = kwds.get("freq", "D")
+
         miss = (checks.missing_any(da, freq) for da in args)
         return reduce(np.logical_or, miss)
 


Test patch:
diff --git a/tests/test_indices.py b/tests/test_indices.py
index 5573b9ce..7df3a3f0 100644
--- a/tests/test_indices.py
+++ b/tests/test_indices.py
@@ -748,6 +748,16 @@ class TestTGXN90p:
         assert out[5] == 25
 
 
+class TestTas:
+    def test_tas(self, tasmin_series, tasmax_series, tas_series):
+        tas = tas_series(np.ones(10))
+        tasmin = tasmin_series(np.zeros(10))
+        tasmax = tasmax_series(np.ones(10) * 2)
+
+        tas_xc = xci.tas(tasmin, tasmax)
+        xr.testing.assert_equal(tas, tas_xc)
+
+
 class TestTxMin:
     @staticmethod
     def time_series(values):
diff --git a/tests/test_temperature.py b/tests/test_temperature.py
index 4cfce560..776a6c77 100644
--- a/tests/test_temperature.py
+++ b/tests/test_temperature.py
@@ -175,13 +175,16 @@ class TestETR:
 
 
 class TestTmean:
-    nc_file = os.path.join(
-        TESTS_DATA, "NRCANdaily", "nrcan_canada_daily_tasmax_1990.nc"
+    nc_files = (
+        os.path.join(TESTS_DATA, "NRCANdaily", "nrcan_canada_daily_tasmax_1990.nc"),
+        os.path.join(TESTS_DATA, "NRCANdaily", "nrcan_canada_daily_tasmin_1990.nc"),
     )
 
     def test_Tmean_3d_data(self):
-        tas = xr.open_dataset(self.nc_file).tasmax
-        tas_C = xr.open_dataset(self.nc_file).tasmax
+        ds_tmax = xr.open_dataset(self.nc_files[0])
+        ds_tmin = xr.open_dataset(self.nc_files[1])
+        tas = atmos.tg(ds_tmin.tasmin, ds_tmax.tasmax)
+        tas_C = atmos.tg(ds_tmin.tasmin, ds_tmax.tasmax)
         tas_C.values -= K2C
         tas_C.attrs["units"] = "C"
         # put a nan somewhere