Repository: xgcm/xgcm
Pull request #270 (2020-10-12 19:09:40)

Emphasize the link to the docs as actionable step
Emphasize the connection to the docs more clearly in the README


Patch:
diff --git a/README.rst b/README.rst
index 98fce41..7580e8d 100644
--- a/README.rst
+++ b/README.rst
@@ -41,8 +41,11 @@ evolving, open-source, python-based framework for GCM analysis will enhance
 the productivity of the field as a whole, accelerating the rate of discovery in
 climate science. xgcm is part of the Pangeo_ initiative.
 
-For more information, including installation instructions, read the full
-`xgcm documentation`_.
+Getting Started
+---------------
+
+To learn how to install and use xgcm for your dataset, visit the `xgcm documentation`_.
+
 
 .. _Pangeo: http://pangeo-data.github.io
 .. _dask: http://dask.pydata.org
diff --git a/doc/whats-new.rst b/doc/whats-new.rst
index 5021f76..bff7b37 100644
--- a/doc/whats-new.rst
+++ b/doc/whats-new.rst
@@ -22,6 +22,8 @@ New Features
   
 Bug fixes
 ~~~~~~~~~
+- More reliable handling of missing values in :py:meth:`Grid.average`. Missing values between data and metrics do not have to be aligned by the user anymore. (:pull:`259`). By `Julius Busecke <https://github.com/jbusecke>`_.
+
 - Remove outdated `example_notebooks` folder (:pull:`244`, :issue:`243`). By `Nikolay Koldunov <https://github.com/koldunovn>`_ and `Julius Busecke <https://github.com/jbusecke>`_.
 .. _whats-new.0.5.0:
 
diff --git a/xgcm/grid.py b/xgcm/grid.py
index faad5e8..8090b76 100644
--- a/xgcm/grid.py
+++ b/xgcm/grid.py
@@ -861,6 +861,31 @@ class Axis:
                 "`transform` can only be used on axes that are non-periodic. Pass `periodic=False` to `xgcm.Grid`."
             )
 
+        def _target_data_name_handling(target_data):
+            """Handle target_data input without a name"""
+            if target_data.name is None:
+                warnings.warn(
+                    "Input`target_data` has no name, but we need a name for the transformed dimension. The name `TRANSFORMED_DIMENSION` will be used. To avoid this warning, call `.rename` on `target_data` before calling `transform`."
+                )
+                target_data.name = "TRANSFORMED_DIMENSION"
+
+        def _check_other_dims(target_da):
+            # check if other dimensions (excluding ones associated with the transform axis) are the
+            # same between `da` and `target_data`. If not provide instructions how to work around.
+
+            da_other_dims = set(da.dims) - set(self.coords.values())
+            target_da_other_dims = set(target_da.dims) - set(self.coords.values())
+            if not target_da_other_dims.issubset(da_other_dims):
+                raise ValueError(
+                    f"Found additional dimensions [{target_da_other_dims-da_other_dims}]"
+                    "in `target_data` not found in `da`. This could mean that the target "
+                    "array is not on the same position along other axes."
+                    " If the additional dimensions are associated witha staggered axis, "
+                    "use grid.interp() to move values to other grid position. "
+                    "If additional dimensions are not related to the grid (e.g. climate "
+                    "model ensemble members or similar), use xr.broadcast() before using transform."
+                )
+
         def _parse_target(target, target_dim, target_data_dim, target_data):
             """Parse target values into correct xarray naming and set default naming based on input data"""
             # if target_data is not provided, assume the target to be one of the staggered dataset dimensions.
@@ -879,29 +904,20 @@ class Axis:
                         )
             else:
                 # if the target is not provided as xr.Dataarray we take the name of the target_data as new dimension name
+                _target_data_name_handling(target_data)
                 target_dim = target_data.name
                 target = xr.DataArray(
                     target, dims=[target_dim], coords={target_dim: target}
                 )
-            return target, target_dim, target_data
 
-        def _check_target_alignment(da, target_da):
-            # check alignment for all other axes, to avoid broadcasting enourmous arrays
-            # to do this on the axis level we will simply check if target_data has some coordinates
-            # that are not found in da (excluding all coordinates that belong to the current axis)
-            da_other_dims = set(da.dims) - set(self.coords.values())
-            target_da_other_dims = set(target_da.dims) - set(self.coords.values())
-            if not target_da_other_dims.issubset(da_other_dims):
-                raise ValueError(
-                    f"Found additional dimensions [{target_da_other_dims-da_other_dims}] in `target_data` not found in `da`. This could mean that the target array is not on the same position along other axes. Use grid.interp() to align the arrays."
-                )
+            _check_other_dims(target_data)
+            return target, target_dim, target_data
 
         _, dim = self._get_axis_coord(da)
         if method == "linear":
             target, target_dim, target_data = _parse_target(
                 target, target_dim, dim, target_data
             )
-            _check_target_alignment(da, target_data)
             out = linear_interpolation(
                 da,
                 target_data,
@@ -927,8 +943,6 @@ class Axis:
                 target, target_dim, target_data_dim, target_data
             )
 
-            _check_target_alignment(da, target_data)
-
             # check on which coordinate `target_data` is, and interpolate if needed
             if target_data_dim not in target_data.dims:
                 warnings.warn(
@@ -1712,24 +1726,20 @@ class Grid:
         axis : str, list of str
             Name of the axis on which to act
         **kwargs: dict
-            Additional arguments passed to `xarray.DataArray.sum`
-
+            Additional arguments passed to `xarray.DataArray.weighted.mean`
 
         Returns
         -------
         da_i : xarray.DataArray
             The averaged data
         """
-
         weight = self.get_metric(da, axis)
-        weighted = da * weight
-        # TODO: We should integrate xarray.weighted once available.
+        weighted = da.weighted(weight)
 
         # get dimension(s) corresponding
         # to `da` and `axis` input
         dim = self._get_dims_from_axis(da, axis)
-        # do we need to pass kwargs?
-        return weighted.sum(dim, **kwargs) / weight.sum(dim, **kwargs)
+        return weighted.mean(dim, **kwargs)
 
     def transform(self, da, axis, target, **kwargs):
         """Convert an array of data to new 1D-coordinates along `axis`.


Test patch:
diff --git a/xgcm/test/test_metric_ops_multi_axis.py b/xgcm/test/test_metric_ops_multi_axis.py
index 5309990..074fca4 100644
--- a/xgcm/test/test_metric_ops_multi_axis.py
+++ b/xgcm/test/test_metric_ops_multi_axis.py
@@ -1,5 +1,7 @@
 from __future__ import print_function
 import pytest
+import xarray as xr
+from xarray.testing import assert_allclose
 
 from xgcm.grid import Grid
 from xgcm.test.datasets import datasets_grid_metric
@@ -50,12 +52,12 @@ class TestParametrized:
             expected = _expected_result(
                 ds.tracer, ds[metric_name], grid, dim, axis, funcname, **kwargs
             )
-            assert new.equals(expected)
+            assert_allclose(new, expected)
 
             # test with tuple input if list is provided
             if isinstance(axis, list):
                 new = func(ds.tracer, tuple(axis), **kwargs)
-                assert new.equals(expected)
+                assert_allclose(new, expected)
 
         # test u position
         for axis, metric_name, dim in zip(
@@ -67,7 +69,7 @@ class TestParametrized:
             expected = _expected_result(
                 ds.u, ds[metric_name], grid, dim, axis, funcname, **kwargs
             )
-            assert new.equals(expected)
+            assert_allclose(new, expected)
 
         # test v position
         for axis, metric_name, dim in zip(
@@ -79,7 +81,7 @@ class TestParametrized:
             expected = _expected_result(
                 ds.v, ds[metric_name], grid, dim, axis, funcname, **kwargs
             )
-            assert new.equals(expected)
+            assert_allclose(new, expected)
 
     def test_cgrid(self, funcname, boundary, periodic):
         ds, coords, metrics = datasets_grid_metric("C")
@@ -104,11 +106,11 @@ class TestParametrized:
             expected = _expected_result(
                 ds.tracer, ds[metric_name], grid, dim, axis, funcname, **kwargs
             )
-            assert new.equals(expected)
+            assert_allclose(new, expected)
             # test with tuple input if list is provided
             if isinstance(axis, list):
                 new = func(ds.tracer, tuple(axis), **kwargs)
-                assert new.equals(expected)
+                assert_allclose(new, expected)
 
         # test u positon
         for axis, metric_name, dim in zip(
@@ -120,7 +122,7 @@ class TestParametrized:
             expected = _expected_result(
                 ds.u, ds[metric_name], grid, dim, axis, funcname, **kwargs
             )
-            assert new.equals(expected)
+            assert_allclose(new, expected)
 
         # test v positon
         for axis, metric_name, dim in zip(
@@ -132,7 +134,7 @@ class TestParametrized:
             expected = _expected_result(
                 ds.v, ds[metric_name], grid, dim, axis, funcname, **kwargs
             )
-            assert new.equals(expected)
+            assert_allclose(new, expected)
 
     @pytest.mark.parametrize("axis", ["X", "Y", "Z"])
     def test_missingaxis(self, axis, funcname, periodic, boundary):
diff --git a/xgcm/test/test_metrics_ops_single_axis.py b/xgcm/test/test_metrics_ops_single_axis.py
index 6d1ffd9..f74b668 100644
--- a/xgcm/test/test_metrics_ops_single_axis.py
+++ b/xgcm/test/test_metrics_ops_single_axis.py
@@ -108,6 +108,27 @@ def test_boundary_global_input(funcname, boundary, fill_value):
     xr.testing.assert_allclose(global_result, manual_result)
 
 
+def test_average_unmatched_missing():
+    # Tests the behavior of grid.average on an array which has missing values, not present in the metric
+    x = np.arange(10)
+    data = xr.DataArray(np.ones(10), dims="x", coords={"x": x})
+    weights = data * 30
+    ds = xr.Dataset({"data": data})
+    ds = ds.assign_coords(weights=weights)
+    # create an xgcm grid
+    grid = Grid(ds, coords={"X": {"center": "x"}}, metrics={"X": ["weights"]})
+
+    # average the unmasked array
+    expected = grid.average(ds.data, "X")
+
+    # now lets introduce a missing value in the data
+    ds.data[6:8] = np.nan
+
+    # assert that the result for both the full and the masked array is equal,
+    # since both only have ones in them.
+    xr.testing.assert_allclose(expected, grid.average(ds.data, "X"))
+
+
 def test_derivative_uniform_grid():
     # this is a uniform grid
     # a non-uniform grid would provide a more rigorous test
diff --git a/xgcm/test/test_transform.py b/xgcm/test/test_transform.py
index 0a006cf..2013190 100644
--- a/xgcm/test/test_transform.py
+++ b/xgcm/test/test_transform.py
@@ -509,9 +509,6 @@ def construct_test_source_data(case_param_dict):
     )
 
 
-# TODO: I am not sure how we would handle periodic axes atm. Should we raise an error?
-
-
 @pytest.fixture(
     # scope="module",
     params=list(cases.keys()),
@@ -795,7 +792,8 @@ def test_conservative_interp_warn():
 
 
 @pytest.mark.skipif(numba is None, reason="numba required")
-def test_grid_transform_noname(multidim_cases):
+def test_grid_transform_noname_data(multidim_cases):
+    """ Check handling of a `da` input without name"""
     source, grid_kwargs, target, transform_kwargs, expected, error_flag = multidim_cases
 
     axis = list(grid_kwargs["coords"].keys())[0]
@@ -810,6 +808,36 @@ def test_grid_transform_noname(multidim_cases):
     assert transformed.name is None
 
 
+@pytest.mark.skipif(numba is None, reason="numba required")
+def test_grid_transform_noname_targetdata():
+    """ Check handling of a `target_data` input without name"""
+    (
+        source,
+        grid_kwargs,
+        target,
+        transform_kwargs,
+        expected,
+        error_flag,
+    ) = construct_test_source_data(cases["linear_depth_dens"])
+
+    axis = list(grid_kwargs["coords"].keys())[0]
+
+    grid = Grid(source, periodic=False, **grid_kwargs)
+
+    source_da = source.data
+    target_data = transform_kwargs.pop("target_data")
+    target_data.name = None
+    # the name of target_data is only used if `target` is provided as numpy array
+    target = target.data
+
+    # the high level routines should be able to deal with all cases (no error flag exception like in the mid level)
+    with pytest.warns(UserWarning):
+        transformed = grid.transform(
+            source_da, axis, target, target_data=target_data, **transform_kwargs
+        )
+    "TRANSFORMED_DIMENSION" in transformed.dims
+
+
 @pytest.mark.skipif(numba is None, reason="numba required")
 def test_transform_error_periodic(multidim_cases):
     source, grid_kwargs, target, transform_kwargs, expected, error_flag = multidim_cases
@@ -922,7 +950,7 @@ def test_grid_transform_multidim(request, client, multidim_cases):
 
 
 @pytest.mark.skipif(numba is None, reason="numba required")
-def test_grid_transform_multidim_alignment_error(request, multidim_cases):
+def test_grid_transform_multidim_other_dims_error(request, multidim_cases):
     # broadcast the 1d column agains some other dims and make sure that the 1d results are still valid
     source, grid_kwargs, target, transform_kwargs, expected, error_flag = multidim_cases