From 59ff68885b3576815843393f047d8eea619958d3 Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Tue, 4 Sep 2018 09:05:39 +0900 Subject: [PATCH 01/24] Added xr.gradient, DataArray.gradient, Dataset.gradient --- xarray/__init__.py | 2 +- xarray/core/computation.py | 62 +++++++++++++++++++++++++++++++++++ xarray/core/dataarray.py | 23 +++++++++++++ xarray/core/dataset.py | 43 ++++++++++++++++++++++-- xarray/core/duck_array_ops.py | 1 + 5 files changed, 127 insertions(+), 4 deletions(-) diff --git a/xarray/__init__.py b/xarray/__init__.py index 7cc7811b783..968b2610350 100644 --- a/xarray/__init__.py +++ b/xarray/__init__.py @@ -10,7 +10,7 @@ from .core.alignment import align, broadcast, broadcast_arrays from .core.common import full_like, zeros_like, ones_like from .core.combine import concat, auto_combine -from .core.computation import apply_ufunc, where, dot +from .core.computation import apply_ufunc, dot, gradient, where from .core.extensions import (register_dataarray_accessor, register_dataset_accessor) from .core.variable import as_variable, Variable, IndexVariable, Coordinate diff --git a/xarray/core/computation.py b/xarray/core/computation.py index bdba72cb48a..cd7331eb511 100644 --- a/xarray/core/computation.py +++ b/xarray/core/computation.py @@ -1137,3 +1137,65 @@ def where(cond, x, y): join='exact', dataset_join='exact', dask='allowed') + + +def _gradient_once(variable, coord, edge_order): + """ Compute the gradient along 1 dimension for Variable. """ + from .variable import Variable + + result_array = duck_array_ops.gradient( + variable.data, coord.data, edge_order=edge_order, + axis=variable.get_axis_num(variable.dims[0])) + return Variable(variable.dims, result_array) + + +def gradient(dataarray, coords, edge_order=1): + """ Compute the gradient of the dataarray with the second order accurate + central differences. + + Parameters + ---------- + dataarray: xr.DataArray + Target array + coords: str or sequence of strs + The coordinates along which the gradient is to be computed. + edge_order: 1 or 2. Default 1 + N-th order accurate differences at the boundaries. + + Returns + ------- + gradient: DataArray or sequence of DataArrays + + See also + -------- + numpy.gradient: corresponding numpy function + """ + from .dataarray import DataArray + + if not isinstance(dataarray, DataArray): + raise TypeError('Only xr.DataArray is supported.' + 'Given {}.'.format([type(arr) for arr in arrays])) + + return_sequence = True + if not isinstance(coords, (tuple, list)): + coords = (coords, ) + return_sequence = False + + result = [] + for coord in coords: + if coord not in dataarray.coords: + raise ValueError('Coordiante {} does not exist.'.format(coord)) + + coord_var = dataarray[coord].variable + if coord_var != 1: + raise ValueError( + 'Only 1d-coordinate is supported. {} is {} ' + 'dimensional.'.format(coord, dataarray[coord].ndim)) + + result.append(DataArray( + _gradient_once(dataarray.variable, coord_var, edge_order), + dims=dataarray.dims, coords=dataarray.coords)) + + if return_sequence: + return tuple(result) + return result[0] diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py index 373a6a4cc9e..3f0fef43254 100644 --- a/xarray/core/dataarray.py +++ b/xarray/core/dataarray.py @@ -2267,6 +2267,29 @@ def rank(self, dim, pct=False, keep_attrs=False): ds = self._to_temp_dataset().rank(dim, pct=pct, keep_attrs=keep_attrs) return self._from_temp_dataset(ds) + def gradient(self, coord, edge_order=1): + """ Compute the gradient with the second order accurate central + differences. + + Parameters + ---------- + coords: str + The coordinate along which the gradient is to be computed. + edge_order: 1 or 2. Default 1 + N-th order accurate differences at the boundaries. + + Returns + ------- + gradient: DataArray + + See also + -------- + numpy.gradient: corresponding numpy function + xr.gradient: more numpy-like function for xarray object. + """ + ds = self._to_temp_dataset().gradient(coord, edge_order) + return self._from_temp_dataset(ds) + # priority most be higher than Variable to properly work with binary ufuncs ops.inject_all_ops_and_reduce_methods(DataArray, priority=60) diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py index 06b258ae261..5278d859faf 100644 --- a/xarray/core/dataset.py +++ b/xarray/core/dataset.py @@ -13,8 +13,8 @@ import xarray as xr from . import ( - alignment, duck_array_ops, formatting, groupby, indexing, ops, resample, - rolling, utils) + alignment, computation, duck_array_ops, formatting, groupby, indexing, ops, + resample, rolling, utils) from .. import conventions from .alignment import align from .common import ( @@ -3645,6 +3645,43 @@ def rank(self, dim, pct=False, keep_attrs=False): attrs = self.attrs if keep_attrs else None return self._replace_vars_and_dims(variables, coord_names, attrs=attrs) + def gradient(self, coord, edge_order=1): + """ Compute the gradient with the second order accurate central + differences. + + Parameters + ---------- + coords: str + The coordinate along which the gradient is to be computed. + edge_order: 1 or 2. Default 1 + N-th order accurate differences at the boundaries. + + Returns + ------- + gradient: Dataset + + See also + -------- + numpy.gradient: corresponding numpy function + xr.gradient: more numpy-like function for xarray object. + """ + if coord not in self.variables: + raise ValueError('Coordinate {} does not exist.'.format(coord)) + + coord_var = self.variables[coord] + if coord_var.ndims: + raise ValueError('Coordinate {} must be 1 dimensional but {} is {}' + ' dimensional'.format(coord, coord_var.ndims)) + + dim = coord_var.dims[0] + variables = OrderedDict() + for k, v in self.variables.items(): + if k in self.data_vars and dim in v.dims: + variables[k] = computation._gradient_once(v, coord_var) + else: + variable[k] = v + return self._replace_vars_and_dims(variables) + @property def real(self): return self._unary_op(lambda x: x.real, keep_attrs=True)(self) @@ -3658,7 +3695,7 @@ def filter_by_attrs(self, **kwargs): Can pass in ``key=value`` or ``key=callable``. A Dataset is returned containing only the variables for which all the filter tests pass. - These tests are either ``key=value`` for which the attribute ``key`` + These tests are either ``key=value`` for which the attribute ``key`` has the exact value ``value`` or the callable passed into ``key=callable`` returns True. The callable will be passed a single value, either the value of the attribute ``key`` or ``None`` if the diff --git a/xarray/core/duck_array_ops.py b/xarray/core/duck_array_ops.py index 17eb310f8db..3520ae02770 100644 --- a/xarray/core/duck_array_ops.py +++ b/xarray/core/duck_array_ops.py @@ -92,6 +92,7 @@ def isnull(data): tensordot = _dask_or_eager_func('tensordot', array_args=slice(2)) einsum = _dask_or_eager_func('einsum', array_args=slice(1, None), requires_dask='0.17.3') +gradient = _dask_or_eager_func('gradient') masked_invalid = _dask_or_eager_func( 'masked_invalid', eager_module=np.ma, From 0adfc68f374f4222e743566477fab0e63cede579 Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Tue, 4 Sep 2018 09:56:40 +0900 Subject: [PATCH 02/24] Working with np.backend --- xarray/core/computation.py | 12 +++++---- xarray/core/dataset.py | 11 +++++---- xarray/tests/test_computation.py | 42 ++++++++++++++++++++++++++++++++ 3 files changed, 55 insertions(+), 10 deletions(-) diff --git a/xarray/core/computation.py b/xarray/core/computation.py index cd7331eb511..a40f9bf87a5 100644 --- a/xarray/core/computation.py +++ b/xarray/core/computation.py @@ -1140,12 +1140,14 @@ def where(cond, x, y): def _gradient_once(variable, coord, edge_order): - """ Compute the gradient along 1 dimension for Variable. """ + """ Compute the gradient along 1 dimension for Variable. + variable, coord: Variable + """ from .variable import Variable result_array = duck_array_ops.gradient( variable.data, coord.data, edge_order=edge_order, - axis=variable.get_axis_num(variable.dims[0])) + axis=variable.get_axis_num(coord.dims[0])) return Variable(variable.dims, result_array) @@ -1183,18 +1185,18 @@ def gradient(dataarray, coords, edge_order=1): result = [] for coord in coords: - if coord not in dataarray.coords: + if coord not in dataarray.coords and coord not in dataarray.dims: raise ValueError('Coordiante {} does not exist.'.format(coord)) coord_var = dataarray[coord].variable - if coord_var != 1: + if coord_var.ndim != 1: raise ValueError( 'Only 1d-coordinate is supported. {} is {} ' 'dimensional.'.format(coord, dataarray[coord].ndim)) result.append(DataArray( _gradient_once(dataarray.variable, coord_var, edge_order), - dims=dataarray.dims, coords=dataarray.coords)) + coords=dataarray.coords)) if return_sequence: return tuple(result) diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py index 5278d859faf..f49a90fbd2e 100644 --- a/xarray/core/dataset.py +++ b/xarray/core/dataset.py @@ -3665,11 +3665,11 @@ def gradient(self, coord, edge_order=1): numpy.gradient: corresponding numpy function xr.gradient: more numpy-like function for xarray object. """ - if coord not in self.variables: + if coord not in self.variables and coord not in self.dims: raise ValueError('Coordinate {} does not exist.'.format(coord)) - coord_var = self.variables[coord] - if coord_var.ndims: + coord_var = self[coord].variable + if coord_var.ndim != 1: raise ValueError('Coordinate {} must be 1 dimensional but {} is {}' ' dimensional'.format(coord, coord_var.ndims)) @@ -3677,9 +3677,10 @@ def gradient(self, coord, edge_order=1): variables = OrderedDict() for k, v in self.variables.items(): if k in self.data_vars and dim in v.dims: - variables[k] = computation._gradient_once(v, coord_var) + variables[k] = computation._gradient_once( + v, coord_var, edge_order) else: - variable[k] = v + variables[k] = v return self._replace_vars_and_dims(variables) @property diff --git a/xarray/tests/test_computation.py b/xarray/tests/test_computation.py index ca8e4e59737..139ed8c2a89 100644 --- a/xarray/tests/test_computation.py +++ b/xarray/tests/test_computation.py @@ -14,6 +14,7 @@ _UFuncSignature, apply_ufunc, broadcast_compat_data, collect_dict_values, join_dict_keys, ordered_set_intersection, ordered_set_union, result_name, unified_dim_sizes) +from xarray.testing import assert_equal from . import raises_regex, requires_dask, has_dask @@ -976,3 +977,44 @@ def test_where(): actual = xr.where(cond, 1, 0) expected = xr.DataArray([1, 0], dims='x') assert_identical(expected, actual) + + +@pytest.mark.parametrize('dask', [True, False]) +@pytest.mark.parametrize('edge_order', [1, 2]) +def test_gradient(dask, edge_order): + rs = np.random.RandomState(42) + da = xr.DataArray(rs.randn(4, 6), dims=['x', 'y'], + coords={'x': [0.4, 0.6, 0.7, 0.75], 'z': 3, + 'x2d': (('x', 'y'), rs.randn(4, 6))}) + if dask and has_dask: + da = da.chunk({'x': 2}) + + ds = xr.Dataset({'var': da}) + + # along x + actual = xr.gradient(da, 'x', edge_order) + expected_x = xr.DataArray( + np.gradient(da, da['x'], axis=0, edge_order=edge_order), + dims=da.dims, coords=da.coords) + assert_equal(expected_x, actual) + assert_equal(actual, ds.gradient('x', edge_order=edge_order)['var']) + assert_equal(ds['var'].gradient('x', edge_order=edge_order), + ds.gradient('x', edge_order=edge_order)['var']) + + # along y + actual = xr.gradient(da, 'y', edge_order) + expected_y = xr.DataArray( + np.gradient(da, da['y'], axis=1, edge_order=edge_order), + dims=da.dims, coords=da.coords) + assert_equal(expected_y, actual) + assert_equal(actual, ds.gradient('y', edge_order=edge_order)['var']) + assert_equal(ds['var'].gradient('y', edge_order=edge_order), + ds.gradient('y', edge_order=edge_order)['var']) + + # along x and y + actual = xr.gradient(da, ('x', 'y'), edge_order) + assert_equal(expected_x, actual[0]) + assert_equal(expected_y, actual[1]) + + with pytest.raises(ValueError): + xr.gradient(da, ('x2d')) From d665b3c636331b097a1a5b205ee42f5e3415185d Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Tue, 4 Sep 2018 14:34:43 +0900 Subject: [PATCH 03/24] test is not passing --- xarray/core/dask_array_ops.py | 46 +++++++++++++++++++++++++++++ xarray/core/duck_array_ops.py | 11 +++++-- xarray/tests/test_computation.py | 8 ++--- xarray/tests/test_duck_array_ops.py | 23 +++++++++++++-- 4 files changed, 79 insertions(+), 9 deletions(-) diff --git a/xarray/core/dask_array_ops.py b/xarray/core/dask_array_ops.py index 423a65aa3c2..4d09df0e94e 100644 --- a/xarray/core/dask_array_ops.py +++ b/xarray/core/dask_array_ops.py @@ -104,3 +104,49 @@ def func(x, window, axis=-1): index = (slice(None),) * axis + (slice(drop_size, drop_size + orig_shape[axis]), ) return out[index] + + +def gradient(a, coord, axis, edge_order): + """ Apply gradient along one dimension """ + if axis < 0: + axis = a.ndim + axis + + for c in a.chunks[axis]: + if c < edge_order + 1: + raise ValueError('Chunk size must be larger than edge_order + 1. ' + 'Given {}. Rechunk to proceed.'.format(c)) + + depth = {d: 1 if d == axis else 0 for d in range(a.ndim)} + # temporary pad zero at the boundary + boundary={d: 0 for d in range(a.ndim)} + ag = overlap(a, depth=depth, boundary=boundary) + + n_chunk = len(a.chunks[axis]) + array_loc_stop = np.cumsum(np.array(a.chunks[axis])) + 2 + array_loc_start = array_loc_stop - np.array(a.chunks[axis]) - 2 + # extrapolate the edge point temporary + coord_pad = np.concatenate([2 * coord[: 1] - coord[1: 2], coord, + 2 * coord[-1:] - coord[-2: -1]]) + + def func(x, block_id): + x = np.asarray(x) + block_loc = block_id[axis] + c = coord_pad[array_loc_start[block_loc]: array_loc_stop[block_loc]] + grad = np.gradient(x, c, axis=axis, edge_order=edge_order) + + # overwrite the boundary value + if block_loc == 0: + idx = (slice(None), ) * axis + (slice(1, edge_order + 2), ) + c1 = c[1: edge_order + 2] + g1 = np.gradient(x[idx], c1, axis=axis, edge_order=edge_order) + grad[(slice(None), ) * axis + (1, )] = g1[ + (slice(None), ) * axis + (1, )] + if block_loc == n_chunk: + idx = (slice(None), ) * axis + (slice(-edge_order - 2, -1), ) + c1 = c[-edge_order - 2: -1] + g1 = np.gradient(x[idx], c1, axis=axis, edge_order=edge_order) + grad[(slice(None), ) * axis + (-2, )] = g1[ + (slice(None), ) * axis + (-2, )] + return grad + + return trim_internal(ag.map_blocks(func, dtype=a.dtype), depth) diff --git a/xarray/core/duck_array_ops.py b/xarray/core/duck_array_ops.py index 3520ae02770..bfeaac71590 100644 --- a/xarray/core/duck_array_ops.py +++ b/xarray/core/duck_array_ops.py @@ -19,7 +19,7 @@ try: import dask.array as dask_array - from . import dask_array_compat + from . import dask_array_compat, dask_array_ops except ImportError: dask_array = None dask_array_compat = None @@ -92,7 +92,14 @@ def isnull(data): tensordot = _dask_or_eager_func('tensordot', array_args=slice(2)) einsum = _dask_or_eager_func('einsum', array_args=slice(1, None), requires_dask='0.17.3') -gradient = _dask_or_eager_func('gradient') + + +def gradient(x, coord, axis, edge_order): + if isinstance(x, dask_array.Array): + return dask_array_ops.gradient( + x, coord, axis=axis, edge_order=edge_order) + return np.gradient(x, coord, axis=axis, edge_order=edge_order) + masked_invalid = _dask_or_eager_func( 'masked_invalid', eager_module=np.ma, diff --git a/xarray/tests/test_computation.py b/xarray/tests/test_computation.py index 139ed8c2a89..495707defba 100644 --- a/xarray/tests/test_computation.py +++ b/xarray/tests/test_computation.py @@ -983,11 +983,11 @@ def test_where(): @pytest.mark.parametrize('edge_order', [1, 2]) def test_gradient(dask, edge_order): rs = np.random.RandomState(42) - da = xr.DataArray(rs.randn(4, 6), dims=['x', 'y'], - coords={'x': [0.4, 0.6, 0.7, 0.75], 'z': 3, - 'x2d': (('x', 'y'), rs.randn(4, 6))}) + da = xr.DataArray(rs.randn(8, 6), dims=['x', 'y'], + coords={'x': [0.2, 0.35, 0.4, 0.6, 0.7, 0.75, 0.76, 0.8], + 'z': 3, 'x2d': (('x', 'y'), rs.randn(8, 6))}) if dask and has_dask: - da = da.chunk({'x': 2}) + da = da.chunk({'x': 4}) ds = xr.Dataset({'var': da}) diff --git a/xarray/tests/test_duck_array_ops.py b/xarray/tests/test_duck_array_ops.py index 3f32fc49fd2..224e00c5f70 100644 --- a/xarray/tests/test_duck_array_ops.py +++ b/xarray/tests/test_duck_array_ops.py @@ -12,8 +12,8 @@ from xarray import DataArray, Dataset, concat from xarray.core import duck_array_ops, dtypes from xarray.core.duck_array_ops import ( - array_notnull_equiv, concatenate, count, first, last, mean, rolling_window, - stack, where) + array_notnull_equiv, concatenate, count, first, gradient, last, mean, + rolling_window, stack, where) from xarray.core.pycompat import dask_array_type from xarray.testing import assert_allclose, assert_equal @@ -258,7 +258,7 @@ def assert_dask_array(da, dask): # TODO test cumsum, cumprod @pytest.mark.parametrize('skipna', [False, True]) @pytest.mark.parametrize('aggdim', [None, 'x']) -def test_reduce(dim_num, dtype, dask, func, skipna, aggdim): +def _test_reduce(dim_num, dtype, dask, func, skipna, aggdim): if aggdim == 'y' and dim_num < 2: pytest.skip('dim not in this test') @@ -417,6 +417,23 @@ def test_dask_rolling(axis, window, center): fill_value=np.nan) +@pytest.mark.skipif(not has_dask, reason='This is for dask.') +@pytest.mark.parametrize('axis', [0, -1, 1]) +@pytest.mark.parametrize('edge_order', [1, 2]) +def test_dask_gradient(axis, edge_order): + import dask.array as da + + array = np.array(np.random.randn(100, 5, 40)) + x = np.exp(np.linspace(0, 1, array.shape[axis])) + + darray = da.from_array(array, chunks=[(6, 30, 30, 20, 14), 5, 8]) + expected = gradient(array, x, axis=axis, edge_order=edge_order) + actual = gradient(darray, x, axis=axis, edge_order=edge_order) + + assert isinstance(actual, da.Array) + assert_array_equal(actual, expected) + + @pytest.mark.parametrize('dim_num', [1, 2]) @pytest.mark.parametrize('dtype', [float, int, np.float32, np.bool_]) @pytest.mark.parametrize('dask', [False, True]) From e0fa5fdb376b790aef6e4dc99ad7f0c38555e0b0 Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Tue, 4 Sep 2018 16:59:03 +0900 Subject: [PATCH 04/24] Docs --- doc/api.rst | 3 ++ doc/whats-new.rst | 7 ++++- xarray/core/computation.py | 46 +++++++++++++++++++++++++++-- xarray/core/dask_array_ops.py | 6 ++-- xarray/core/dataarray.py | 27 ++++++++++++++++- xarray/core/dataset.py | 2 +- xarray/tests/test_duck_array_ops.py | 2 +- 7 files changed, 84 insertions(+), 9 deletions(-) diff --git a/doc/api.rst b/doc/api.rst index 927c0aa072c..67e396f656f 100644 --- a/doc/api.rst +++ b/doc/api.rst @@ -25,6 +25,7 @@ Top-level functions zeros_like ones_like dot + gradient Dataset ======= @@ -150,6 +151,7 @@ Computation Dataset.resample Dataset.diff Dataset.quantile + Dataset.gradient **Aggregation**: :py:attr:`~Dataset.all` @@ -317,6 +319,7 @@ Computation DataArray.diff DataArray.dot DataArray.quantile + DataArray.gradient **Aggregation**: :py:attr:`~DataArray.all` diff --git a/doc/whats-new.rst b/doc/whats-new.rst index dd9eb9e48fe..ea7176ded5e 100644 --- a/doc/whats-new.rst +++ b/doc/whats-new.rst @@ -36,6 +36,11 @@ Documentation Enhancements ~~~~~~~~~~~~ +- :py:func:`~xarray.gradient`, :py:meth:`~xarray.DataArray.gradient`, and + :py:meth:`~xarray.Dataset.gradient` are newly added. + (:issue:`1332`) + By `Keisuke Fujii `_. + - min_count option is newly supported in :py:meth:`~xarray.DataArray.sum`, :py:meth:`~xarray.DataArray.prod` and :py:meth:`~xarray.Dataset.sum`, and :py:meth:`~xarray.Dataset.prod`. @@ -84,7 +89,7 @@ Bug fixes - Tests can be run in parallel with pytest-xdist By `Tony Tung `_. -- Follow up the renamings in dask; from dask.ghost to dask.overlap +- Follow up the renamings in dask; from dask.ghost to dask.overlap By `Keisuke Fujii `_. - Now raises a ValueError when there is a conflict between dimension names and diff --git a/xarray/core/computation.py b/xarray/core/computation.py index a40f9bf87a5..7ba67f6b647 100644 --- a/xarray/core/computation.py +++ b/xarray/core/computation.py @@ -1160,17 +1160,59 @@ def gradient(dataarray, coords, edge_order=1): dataarray: xr.DataArray Target array coords: str or sequence of strs - The coordinates along which the gradient is to be computed. + The coordinate to be used to compute the gradient. edge_order: 1 or 2. Default 1 N-th order accurate differences at the boundaries. Returns ------- - gradient: DataArray or sequence of DataArrays + gradient: DataArray or a sequence of DataArrays See also -------- numpy.gradient: corresponding numpy function + + Examples + -------- + + >>> da = xr.DataArray(np.arange(12).reshape(4, 3), dims=['x', 'y'], + ... coords={'x': [0, 0.1, 1.1, 1.2]}) + >>> da + + array([[ 0, 1, 2], + [ 3, 4, 5], + [ 6, 7, 8], + [ 9, 10, 11]]) + Coordinates: + * x (x) float64 0.0 0.1 1.1 1.2 + Dimensions without coordinates: y + >>> + >>> xr.gradient(da, 'x') + + array([[30. , 30. , 30. ], + [27.545455, 27.545455, 27.545455], + [27.545455, 27.545455, 27.545455], + [30. , 30. , 30. ]]) + Coordinates: + * x (x) float64 0.0 0.1 1.1 1.2 + Dimensions without coordinates: y + >>> + >>> xr.gradient(da, ('x', 'y')) + ( + array([[30. , 30. , 30. ], + [27.545455, 27.545455, 27.545455], + [27.545455, 27.545455, 27.545455], + [30. , 30. , 30. ]]) + Coordinates: + * x (x) float64 0.0 0.1 1.1 1.2 + Dimensions without coordinates: y, + array([[1., 1., 1.], + [1., 1., 1.], + [1., 1., 1.], + [1., 1., 1.]]) + Coordinates: + * x (x) float64 0.0 0.1 1.1 1.2 + Dimensions without coordinates: y) """ from .dataarray import DataArray diff --git a/xarray/core/dask_array_ops.py b/xarray/core/dask_array_ops.py index 4d09df0e94e..df0d68b7ef3 100644 --- a/xarray/core/dask_array_ops.py +++ b/xarray/core/dask_array_ops.py @@ -140,13 +140,13 @@ def func(x, block_id): c1 = c[1: edge_order + 2] g1 = np.gradient(x[idx], c1, axis=axis, edge_order=edge_order) grad[(slice(None), ) * axis + (1, )] = g1[ - (slice(None), ) * axis + (1, )] - if block_loc == n_chunk: + (slice(None), ) * axis + (0, )] + if block_loc == n_chunk - 1: idx = (slice(None), ) * axis + (slice(-edge_order - 2, -1), ) c1 = c[-edge_order - 2: -1] g1 = np.gradient(x[idx], c1, axis=axis, edge_order=edge_order) grad[(slice(None), ) * axis + (-2, )] = g1[ - (slice(None), ) * axis + (-2, )] + (slice(None), ) * axis + (-1, )] return grad return trim_internal(ag.map_blocks(func, dtype=a.dtype), depth) diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py index 3f0fef43254..3defb37d2cf 100644 --- a/xarray/core/dataarray.py +++ b/xarray/core/dataarray.py @@ -2274,7 +2274,7 @@ def gradient(self, coord, edge_order=1): Parameters ---------- coords: str - The coordinate along which the gradient is to be computed. + The coordinate to be used to compute the gradient. edge_order: 1 or 2. Default 1 N-th order accurate differences at the boundaries. @@ -2286,6 +2286,31 @@ def gradient(self, coord, edge_order=1): -------- numpy.gradient: corresponding numpy function xr.gradient: more numpy-like function for xarray object. + + Examples + -------- + + >>> da = xr.DataArray(np.arange(12).reshape(4, 3), dims=['x', 'y'], + ... coords={'x': [0, 0.1, 1.1, 1.2]}) + >>> da + + array([[ 0, 1, 2], + [ 3, 4, 5], + [ 6, 7, 8], + [ 9, 10, 11]]) + Coordinates: + * x (x) float64 0.0 0.1 1.1 1.2 + Dimensions without coordinates: y + >>> + >>> da.gradient('x') + + array([[30. , 30. , 30. ], + [27.545455, 27.545455, 27.545455], + [27.545455, 27.545455, 27.545455], + [30. , 30. , 30. ]]) + Coordinates: + * x (x) float64 0.0 0.1 1.1 1.2 + Dimensions without coordinates: y """ ds = self._to_temp_dataset().gradient(coord, edge_order) return self._from_temp_dataset(ds) diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py index f49a90fbd2e..d7e4589360a 100644 --- a/xarray/core/dataset.py +++ b/xarray/core/dataset.py @@ -3652,7 +3652,7 @@ def gradient(self, coord, edge_order=1): Parameters ---------- coords: str - The coordinate along which the gradient is to be computed. + The coordinate to be used to compute the gradient. edge_order: 1 or 2. Default 1 N-th order accurate differences at the boundaries. diff --git a/xarray/tests/test_duck_array_ops.py b/xarray/tests/test_duck_array_ops.py index 224e00c5f70..b9712f60290 100644 --- a/xarray/tests/test_duck_array_ops.py +++ b/xarray/tests/test_duck_array_ops.py @@ -258,7 +258,7 @@ def assert_dask_array(da, dask): # TODO test cumsum, cumprod @pytest.mark.parametrize('skipna', [False, True]) @pytest.mark.parametrize('aggdim', [None, 'x']) -def _test_reduce(dim_num, dtype, dask, func, skipna, aggdim): +def test_reduce(dim_num, dtype, dask, func, skipna, aggdim): if aggdim == 'y' and dim_num < 2: pytest.skip('dim not in this test') From 218e62d7507de6c5864db1486a351b5ce69c2408 Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Tue, 4 Sep 2018 17:06:12 +0900 Subject: [PATCH 05/24] flake8 --- xarray/core/computation.py | 4 ++-- xarray/core/dask_array_ops.py | 6 +++--- xarray/core/duck_array_ops.py | 2 +- xarray/tests/test_computation.py | 2 +- 4 files changed, 7 insertions(+), 7 deletions(-) diff --git a/xarray/core/computation.py b/xarray/core/computation.py index 7ba67f6b647..5a9935d145d 100644 --- a/xarray/core/computation.py +++ b/xarray/core/computation.py @@ -1217,8 +1217,8 @@ def gradient(dataarray, coords, edge_order=1): from .dataarray import DataArray if not isinstance(dataarray, DataArray): - raise TypeError('Only xr.DataArray is supported.' - 'Given {}.'.format([type(arr) for arr in arrays])) + raise TypeError( + 'Only DataArray is supported. Given {}.'.format(type(dataarray))) return_sequence = True if not isinstance(coords, (tuple, list)): diff --git a/xarray/core/dask_array_ops.py b/xarray/core/dask_array_ops.py index df0d68b7ef3..66eff43815e 100644 --- a/xarray/core/dask_array_ops.py +++ b/xarray/core/dask_array_ops.py @@ -118,7 +118,7 @@ def gradient(a, coord, axis, edge_order): depth = {d: 1 if d == axis else 0 for d in range(a.ndim)} # temporary pad zero at the boundary - boundary={d: 0 for d in range(a.ndim)} + boundary = {d: 0 for d in range(a.ndim)} ag = overlap(a, depth=depth, boundary=boundary) n_chunk = len(a.chunks[axis]) @@ -140,13 +140,13 @@ def func(x, block_id): c1 = c[1: edge_order + 2] g1 = np.gradient(x[idx], c1, axis=axis, edge_order=edge_order) grad[(slice(None), ) * axis + (1, )] = g1[ - (slice(None), ) * axis + (0, )] + (slice(None), ) * axis + (0, )] if block_loc == n_chunk - 1: idx = (slice(None), ) * axis + (slice(-edge_order - 2, -1), ) c1 = c[-edge_order - 2: -1] g1 = np.gradient(x[idx], c1, axis=axis, edge_order=edge_order) grad[(slice(None), ) * axis + (-2, )] = g1[ - (slice(None), ) * axis + (-1, )] + (slice(None), ) * axis + (-1, )] return grad return trim_internal(ag.map_blocks(func, dtype=a.dtype), depth) diff --git a/xarray/core/duck_array_ops.py b/xarray/core/duck_array_ops.py index bfeaac71590..5676990b78f 100644 --- a/xarray/core/duck_array_ops.py +++ b/xarray/core/duck_array_ops.py @@ -19,7 +19,7 @@ try: import dask.array as dask_array - from . import dask_array_compat, dask_array_ops + from . import dask_array_compat except ImportError: dask_array = None dask_array_compat = None diff --git a/xarray/tests/test_computation.py b/xarray/tests/test_computation.py index 495707defba..364cdd0be85 100644 --- a/xarray/tests/test_computation.py +++ b/xarray/tests/test_computation.py @@ -999,7 +999,7 @@ def test_gradient(dask, edge_order): assert_equal(expected_x, actual) assert_equal(actual, ds.gradient('x', edge_order=edge_order)['var']) assert_equal(ds['var'].gradient('x', edge_order=edge_order), - ds.gradient('x', edge_order=edge_order)['var']) + ds.gradient('x', edge_order=edge_order)['var']) # along y actual = xr.gradient(da, 'y', edge_order) From 888b924883d514f1e931e7343df3ef6ad22cd0d4 Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Tue, 4 Sep 2018 18:07:31 +0900 Subject: [PATCH 06/24] support environment without dask --- xarray/core/duck_array_ops.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/xarray/core/duck_array_ops.py b/xarray/core/duck_array_ops.py index 5676990b78f..d3a6067f542 100644 --- a/xarray/core/duck_array_ops.py +++ b/xarray/core/duck_array_ops.py @@ -95,7 +95,7 @@ def isnull(data): def gradient(x, coord, axis, edge_order): - if isinstance(x, dask_array.Array): + if isinstance(x, dask_array_type): return dask_array_ops.gradient( x, coord, axis=axis, edge_order=edge_order) return np.gradient(x, coord, axis=axis, edge_order=edge_order) From a0ab4c222858615034b16b93199bd1eaad4c744f Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Tue, 4 Sep 2018 20:31:26 +0900 Subject: [PATCH 07/24] Support numpy < 1.13 --- xarray/core/duck_array_ops.py | 2 +- xarray/core/npcompat.py | 316 ++++++++++++++++++++++++++++++++++ 2 files changed, 317 insertions(+), 1 deletion(-) diff --git a/xarray/core/duck_array_ops.py b/xarray/core/duck_array_ops.py index d3a6067f542..0d44c4f108d 100644 --- a/xarray/core/duck_array_ops.py +++ b/xarray/core/duck_array_ops.py @@ -98,7 +98,7 @@ def gradient(x, coord, axis, edge_order): if isinstance(x, dask_array_type): return dask_array_ops.gradient( x, coord, axis=axis, edge_order=edge_order) - return np.gradient(x, coord, axis=axis, edge_order=edge_order) + return npcompat.gradient(x, coord, axis=axis, edge_order=edge_order) masked_invalid = _dask_or_eager_func( diff --git a/xarray/core/npcompat.py b/xarray/core/npcompat.py index 6d4db063b98..98ddc70f37a 100644 --- a/xarray/core/npcompat.py +++ b/xarray/core/npcompat.py @@ -1,5 +1,6 @@ from __future__ import absolute_import, division, print_function +from distutils.version import LooseVersion import numpy as np try: @@ -97,3 +98,318 @@ def isin(element, test_elements, assume_unique=False, invert=False): element = np.asarray(element) return np.in1d(element, test_elements, assume_unique=assume_unique, invert=invert).reshape(element.shape) + + +if LooseVersion(np.__version__) >= LooseVersion('1.13'): + gradient = np.gradient +else: + import numpy.core.numeric as _nx + + def gradient(f, *varargs, **kwargs): + """ + Return the gradient of an N-dimensional array. + The gradient is computed using second order accurate central differences + in the interior points and either first or second order accurate one-sides + (forward or backwards) differences at the boundaries. + The returned gradient hence has the same shape as the input array. + Parameters + ---------- + f : array_like + An N-dimensional array containing samples of a scalar function. + varargs : list of scalar or array, optional + Spacing between f values. Default unitary spacing for all dimensions. + Spacing can be specified using: + 1. single scalar to specify a sample distance for all dimensions. + 2. N scalars to specify a constant sample distance for each dimension. + i.e. `dx`, `dy`, `dz`, ... + 3. N arrays to specify the coordinates of the values along each + dimension of F. The length of the array must match the size of + the corresponding dimension + 4. Any combination of N scalars/arrays with the meaning of 2. and 3. + If `axis` is given, the number of varargs must equal the number of axes. + Default: 1. + edge_order : {1, 2}, optional + Gradient is calculated using N-th order accurate differences + at the boundaries. Default: 1. + .. versionadded:: 1.9.1 + axis : None or int or tuple of ints, optional + Gradient is calculated only along the given axis or axes + The default (axis = None) is to calculate the gradient for all the axes + of the input array. axis may be negative, in which case it counts from + the last to the first axis. + .. versionadded:: 1.11.0 + Returns + ------- + gradient : ndarray or list of ndarray + A set of ndarrays (or a single ndarray if there is only one dimension) + corresponding to the derivatives of f with respect to each dimension. + Each derivative has the same shape as f. + Examples + -------- + >>> f = np.array([1, 2, 4, 7, 11, 16], dtype=np.float) + >>> np.gradient(f) + array([ 1. , 1.5, 2.5, 3.5, 4.5, 5. ]) + >>> np.gradient(f, 2) + array([ 0.5 , 0.75, 1.25, 1.75, 2.25, 2.5 ]) + Spacing can be also specified with an array that represents the coordinates + of the values F along the dimensions. + For instance a uniform spacing: + >>> x = np.arange(f.size) + >>> np.gradient(f, x) + array([ 1. , 1.5, 2.5, 3.5, 4.5, 5. ]) + Or a non uniform one: + >>> x = np.array([0., 1., 1.5, 3.5, 4., 6.], dtype=np.float) + >>> np.gradient(f, x) + array([ 1. , 3. , 3.5, 6.7, 6.9, 2.5]) + For two dimensional arrays, the return will be two arrays ordered by + axis. In this example the first array stands for the gradient in + rows and the second one in columns direction: + >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=np.float)) + [array([[ 2., 2., -1.], + [ 2., 2., -1.]]), array([[ 1. , 2.5, 4. ], + [ 1. , 1. , 1. ]])] + In this example the spacing is also specified: + uniform for axis=0 and non uniform for axis=1 + >>> dx = 2. + >>> y = [1., 1.5, 3.5] + >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=np.float), dx, y) + [array([[ 1. , 1. , -0.5], + [ 1. , 1. , -0.5]]), array([[ 2. , 2. , 2. ], + [ 2. , 1.7, 0.5]])] + It is possible to specify how boundaries are treated using `edge_order` + >>> x = np.array([0, 1, 2, 3, 4]) + >>> f = x**2 + >>> np.gradient(f, edge_order=1) + array([ 1., 2., 4., 6., 7.]) + >>> np.gradient(f, edge_order=2) + array([-0., 2., 4., 6., 8.]) + The `axis` keyword can be used to specify a subset of axes of which the + gradient is calculated + >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=np.float), axis=0) + array([[ 2., 2., -1.], + [ 2., 2., -1.]]) + Notes + ----- + Assuming that :math:`f\\in C^{3}` (i.e., :math:`f` has at least 3 continuous + derivatives) and let be :math:`h_{*}` a non homogeneous stepsize, the + spacing the finite difference coefficients are computed by minimising + the consistency error :math:`\\eta_{i}`: + .. math:: + \\eta_{i} = f_{i}^{\\left(1\\right)} - + \\left[ \\alpha f\\left(x_{i}\\right) + + \\beta f\\left(x_{i} + h_{d}\\right) + + \\gamma f\\left(x_{i}-h_{s}\\right) + \\right] + By substituting :math:`f(x_{i} + h_{d})` and :math:`f(x_{i} - h_{s})` + with their Taylor series expansion, this translates into solving + the following the linear system: + .. math:: + \\left\\{ + \\begin{array}{r} + \\alpha+\\beta+\\gamma=0 \\\\ + -\\beta h_{d}+\\gamma h_{s}=1 \\\\ + \\beta h_{d}^{2}+\\gamma h_{s}^{2}=0 + \\end{array} + \\right. + The resulting approximation of :math:`f_{i}^{(1)}` is the following: + .. math:: + \\hat f_{i}^{(1)} = + \\frac{ + h_{s}^{2}f\\left(x_{i} + h_{d}\\right) + + \\left(h_{d}^{2} - h_{s}^{2}\\right)f\\left(x_{i}\\right) + - h_{d}^{2}f\\left(x_{i}-h_{s}\\right)} + { h_{s}h_{d}\\left(h_{d} + h_{s}\\right)} + + \\mathcal{O}\\left(\\frac{h_{d}h_{s}^{2} + + h_{s}h_{d}^{2}}{h_{d} + + h_{s}}\\right) + It is worth noting that if :math:`h_{s}=h_{d}` + (i.e., data are evenly spaced) + we find the standard second order approximation: + .. math:: + \\hat f_{i}^{(1)}= + \\frac{f\\left(x_{i+1}\\right) - f\\left(x_{i-1}\\right)}{2h} + + \\mathcal{O}\\left(h^{2}\\right) + With a similar procedure the forward/backward approximations used for + boundaries can be derived. + References + ---------- + .. [1] Quarteroni A., Sacco R., Saleri F. (2007) Numerical Mathematics + (Texts in Applied Mathematics). New York: Springer. + .. [2] Durran D. R. (1999) Numerical Methods for Wave Equations + in Geophysical Fluid Dynamics. New York: Springer. + .. [3] Fornberg B. (1988) Generation of Finite Difference Formulas on + Arbitrarily Spaced Grids, + Mathematics of Computation 51, no. 184 : 699-706. + `PDF `_. + """ + f = np.asanyarray(f) + N = f.ndim # number of dimensions + + axes = kwargs.pop('axis', None) + if axes is None: + axes = tuple(range(N)) + else: + axes = _nx.normalize_axis_tuple(axes, N) + + len_axes = len(axes) + n = len(varargs) + if n == 0: + # no spacing argument - use 1 in all axes + dx = [1.0] * len_axes + elif n == 1 and np.ndim(varargs[0]) == 0: + # single scalar for all axes + dx = varargs * len_axes + elif n == len_axes: + # scalar or 1d array for each axis + dx = list(varargs) + for i, distances in enumerate(dx): + if np.ndim(distances) == 0: + continue + elif np.ndim(distances) != 1: + raise ValueError("distances must be either scalars or 1d") + if len(distances) != f.shape[axes[i]]: + raise ValueError("when 1d, distances must match " + "the length of the corresponding dimension") + diffx = np.diff(distances) + # if distances are constant reduce to the scalar case + # since it brings a consistent speedup + if (diffx == diffx[0]).all(): + diffx = diffx[0] + dx[i] = diffx + else: + raise TypeError("invalid number of arguments") + + edge_order = kwargs.pop('edge_order', 1) + if kwargs: + raise TypeError('"{}" are not valid keyword arguments.'.format( + '", "'.join(kwargs.keys()))) + if edge_order > 2: + raise ValueError("'edge_order' greater than 2 not supported") + + # use central differences on interior and one-sided differences on the + # endpoints. This preserves second order-accuracy over the full domain. + + outvals = [] + + # create slice objects --- initially all are [:, :, ..., :] + slice1 = [slice(None)]*N + slice2 = [slice(None)]*N + slice3 = [slice(None)]*N + slice4 = [slice(None)]*N + + otype = f.dtype.char + if otype not in ['f', 'd', 'F', 'D', 'm', 'M']: + otype = 'd' + + # Difference of datetime64 elements results in timedelta64 + if otype == 'M': + # Need to use the full dtype name because it contains unit information + otype = f.dtype.name.replace('datetime', 'timedelta') + elif otype == 'm': + # Needs to keep the specific units, can't be a general unit + otype = f.dtype + + # Convert datetime64 data into ints. Make dummy variable `y` + # that is a view of ints if the data is datetime64, otherwise + # just set y equal to the array `f`. + if f.dtype.char in ["M", "m"]: + y = f.view('int64') + else: + y = f + + for i, axis in enumerate(axes): + if y.shape[axis] < edge_order + 1: + raise ValueError( + "Shape of array too small to calculate a numerical gradient, " + "at least (edge_order + 1) elements are required.") + # result allocation + out = np.empty_like(y, dtype=otype) + + uniform_spacing = np.ndim(dx[i]) == 0 + + # Numerical differentiation: 2nd order interior + slice1[axis] = slice(1, -1) + slice2[axis] = slice(None, -2) + slice3[axis] = slice(1, -1) + slice4[axis] = slice(2, None) + + if uniform_spacing: + out[slice1] = (f[slice4] - f[slice2]) / (2. * dx[i]) + else: + dx1 = dx[i][0:-1] + dx2 = dx[i][1:] + a = -(dx2)/(dx1 * (dx1 + dx2)) + b = (dx2 - dx1) / (dx1 * dx2) + c = dx1 / (dx2 * (dx1 + dx2)) + # fix the shape for broadcasting + shape = np.ones(N, dtype=int) + shape[axis] = -1 + a.shape = b.shape = c.shape = shape + # 1D equivalent -- out[1:-1] = a * f[:-2] + b * f[1:-1] + c * f[2:] + out[slice1] = a * f[slice2] + b * f[slice3] + c * f[slice4] + + # Numerical differentiation: 1st order edges + if edge_order == 1: + slice1[axis] = 0 + slice2[axis] = 1 + slice3[axis] = 0 + dx_0 = dx[i] if uniform_spacing else dx[i][0] + # 1D equivalent -- out[0] = (y[1] - y[0]) / (x[1] - x[0]) + out[slice1] = (y[slice2] - y[slice3]) / dx_0 + + slice1[axis] = -1 + slice2[axis] = -1 + slice3[axis] = -2 + dx_n = dx[i] if uniform_spacing else dx[i][-1] + # 1D equivalent -- out[-1] = (y[-1] - y[-2]) / (x[-1] - x[-2]) + out[slice1] = (y[slice2] - y[slice3]) / dx_n + + # Numerical differentiation: 2nd order edges + else: + slice1[axis] = 0 + slice2[axis] = 0 + slice3[axis] = 1 + slice4[axis] = 2 + if uniform_spacing: + a = -1.5 / dx[i] + b = 2. / dx[i] + c = -0.5 / dx[i] + else: + dx1 = dx[i][0] + dx2 = dx[i][1] + a = -(2. * dx1 + dx2)/(dx1 * (dx1 + dx2)) + b = (dx1 + dx2) / (dx1 * dx2) + c = - dx1 / (dx2 * (dx1 + dx2)) + # 1D equivalent -- out[0] = a * y[0] + b * y[1] + c * y[2] + out[slice1] = a * y[slice2] + b * y[slice3] + c * y[slice4] + + slice1[axis] = -1 + slice2[axis] = -3 + slice3[axis] = -2 + slice4[axis] = -1 + if uniform_spacing: + a = 0.5 / dx[i] + b = -2. / dx[i] + c = 1.5 / dx[i] + else: + dx1 = dx[i][-2] + dx2 = dx[i][-1] + a = (dx2) / (dx1 * (dx1 + dx2)) + b = - (dx2 + dx1) / (dx1 * dx2) + c = (2. * dx2 + dx1) / (dx2 * (dx1 + dx2)) + # 1D equivalent -- out[-1] = a * f[-3] + b * f[-2] + c * f[-1] + out[slice1] = a * y[slice2] + b * y[slice3] + c * y[slice4] + + outvals.append(out) + + # reset the slice object in this dimension to ":" + slice1[axis] = slice(None) + slice2[axis] = slice(None) + slice3[axis] = slice(None) + slice4[axis] = slice(None) + + if len_axes == 1: + return outvals[0] + else: + return outvals From c58151359a32b343ccb03832e9b0e98c9c6f935d Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Wed, 5 Sep 2018 07:18:08 +0900 Subject: [PATCH 08/24] Support numpy 1.12 --- xarray/core/npcompat.py | 173 +++++----------------------------------- 1 file changed, 21 insertions(+), 152 deletions(-) diff --git a/xarray/core/npcompat.py b/xarray/core/npcompat.py index 98ddc70f37a..22dff44acf8 100644 --- a/xarray/core/npcompat.py +++ b/xarray/core/npcompat.py @@ -103,146 +103,12 @@ def isin(element, test_elements, assume_unique=False, invert=False): if LooseVersion(np.__version__) >= LooseVersion('1.13'): gradient = np.gradient else: - import numpy.core.numeric as _nx + def normalize_axis_tuple(axes, N): + if isinstance(axes, int): + axes = (axes, ) + return tuple([N + a if a < 0 else a for a in axes]) def gradient(f, *varargs, **kwargs): - """ - Return the gradient of an N-dimensional array. - The gradient is computed using second order accurate central differences - in the interior points and either first or second order accurate one-sides - (forward or backwards) differences at the boundaries. - The returned gradient hence has the same shape as the input array. - Parameters - ---------- - f : array_like - An N-dimensional array containing samples of a scalar function. - varargs : list of scalar or array, optional - Spacing between f values. Default unitary spacing for all dimensions. - Spacing can be specified using: - 1. single scalar to specify a sample distance for all dimensions. - 2. N scalars to specify a constant sample distance for each dimension. - i.e. `dx`, `dy`, `dz`, ... - 3. N arrays to specify the coordinates of the values along each - dimension of F. The length of the array must match the size of - the corresponding dimension - 4. Any combination of N scalars/arrays with the meaning of 2. and 3. - If `axis` is given, the number of varargs must equal the number of axes. - Default: 1. - edge_order : {1, 2}, optional - Gradient is calculated using N-th order accurate differences - at the boundaries. Default: 1. - .. versionadded:: 1.9.1 - axis : None or int or tuple of ints, optional - Gradient is calculated only along the given axis or axes - The default (axis = None) is to calculate the gradient for all the axes - of the input array. axis may be negative, in which case it counts from - the last to the first axis. - .. versionadded:: 1.11.0 - Returns - ------- - gradient : ndarray or list of ndarray - A set of ndarrays (or a single ndarray if there is only one dimension) - corresponding to the derivatives of f with respect to each dimension. - Each derivative has the same shape as f. - Examples - -------- - >>> f = np.array([1, 2, 4, 7, 11, 16], dtype=np.float) - >>> np.gradient(f) - array([ 1. , 1.5, 2.5, 3.5, 4.5, 5. ]) - >>> np.gradient(f, 2) - array([ 0.5 , 0.75, 1.25, 1.75, 2.25, 2.5 ]) - Spacing can be also specified with an array that represents the coordinates - of the values F along the dimensions. - For instance a uniform spacing: - >>> x = np.arange(f.size) - >>> np.gradient(f, x) - array([ 1. , 1.5, 2.5, 3.5, 4.5, 5. ]) - Or a non uniform one: - >>> x = np.array([0., 1., 1.5, 3.5, 4., 6.], dtype=np.float) - >>> np.gradient(f, x) - array([ 1. , 3. , 3.5, 6.7, 6.9, 2.5]) - For two dimensional arrays, the return will be two arrays ordered by - axis. In this example the first array stands for the gradient in - rows and the second one in columns direction: - >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=np.float)) - [array([[ 2., 2., -1.], - [ 2., 2., -1.]]), array([[ 1. , 2.5, 4. ], - [ 1. , 1. , 1. ]])] - In this example the spacing is also specified: - uniform for axis=0 and non uniform for axis=1 - >>> dx = 2. - >>> y = [1., 1.5, 3.5] - >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=np.float), dx, y) - [array([[ 1. , 1. , -0.5], - [ 1. , 1. , -0.5]]), array([[ 2. , 2. , 2. ], - [ 2. , 1.7, 0.5]])] - It is possible to specify how boundaries are treated using `edge_order` - >>> x = np.array([0, 1, 2, 3, 4]) - >>> f = x**2 - >>> np.gradient(f, edge_order=1) - array([ 1., 2., 4., 6., 7.]) - >>> np.gradient(f, edge_order=2) - array([-0., 2., 4., 6., 8.]) - The `axis` keyword can be used to specify a subset of axes of which the - gradient is calculated - >>> np.gradient(np.array([[1, 2, 6], [3, 4, 5]], dtype=np.float), axis=0) - array([[ 2., 2., -1.], - [ 2., 2., -1.]]) - Notes - ----- - Assuming that :math:`f\\in C^{3}` (i.e., :math:`f` has at least 3 continuous - derivatives) and let be :math:`h_{*}` a non homogeneous stepsize, the - spacing the finite difference coefficients are computed by minimising - the consistency error :math:`\\eta_{i}`: - .. math:: - \\eta_{i} = f_{i}^{\\left(1\\right)} - - \\left[ \\alpha f\\left(x_{i}\\right) + - \\beta f\\left(x_{i} + h_{d}\\right) + - \\gamma f\\left(x_{i}-h_{s}\\right) - \\right] - By substituting :math:`f(x_{i} + h_{d})` and :math:`f(x_{i} - h_{s})` - with their Taylor series expansion, this translates into solving - the following the linear system: - .. math:: - \\left\\{ - \\begin{array}{r} - \\alpha+\\beta+\\gamma=0 \\\\ - -\\beta h_{d}+\\gamma h_{s}=1 \\\\ - \\beta h_{d}^{2}+\\gamma h_{s}^{2}=0 - \\end{array} - \\right. - The resulting approximation of :math:`f_{i}^{(1)}` is the following: - .. math:: - \\hat f_{i}^{(1)} = - \\frac{ - h_{s}^{2}f\\left(x_{i} + h_{d}\\right) - + \\left(h_{d}^{2} - h_{s}^{2}\\right)f\\left(x_{i}\\right) - - h_{d}^{2}f\\left(x_{i}-h_{s}\\right)} - { h_{s}h_{d}\\left(h_{d} + h_{s}\\right)} - + \\mathcal{O}\\left(\\frac{h_{d}h_{s}^{2} - + h_{s}h_{d}^{2}}{h_{d} - + h_{s}}\\right) - It is worth noting that if :math:`h_{s}=h_{d}` - (i.e., data are evenly spaced) - we find the standard second order approximation: - .. math:: - \\hat f_{i}^{(1)}= - \\frac{f\\left(x_{i+1}\\right) - f\\left(x_{i-1}\\right)}{2h} - + \\mathcal{O}\\left(h^{2}\\right) - With a similar procedure the forward/backward approximations used for - boundaries can be derived. - References - ---------- - .. [1] Quarteroni A., Sacco R., Saleri F. (2007) Numerical Mathematics - (Texts in Applied Mathematics). New York: Springer. - .. [2] Durran D. R. (1999) Numerical Methods for Wave Equations - in Geophysical Fluid Dynamics. New York: Springer. - .. [3] Fornberg B. (1988) Generation of Finite Difference Formulas on - Arbitrarily Spaced Grids, - Mathematics of Computation 51, no. 184 : 699-706. - `PDF `_. - """ f = np.asanyarray(f) N = f.ndim # number of dimensions @@ -250,7 +116,7 @@ def gradient(f, *varargs, **kwargs): if axes is None: axes = tuple(range(N)) else: - axes = _nx.normalize_axis_tuple(axes, N) + axes = normalize_axis_tuple(axes, N) len_axes = len(axes) n = len(varargs) @@ -269,8 +135,8 @@ def gradient(f, *varargs, **kwargs): elif np.ndim(distances) != 1: raise ValueError("distances must be either scalars or 1d") if len(distances) != f.shape[axes[i]]: - raise ValueError("when 1d, distances must match " - "the length of the corresponding dimension") + raise ValueError("when 1d, distances must match the " + "length of the corresponding dimension") diffx = np.diff(distances) # if distances are constant reduce to the scalar case # since it brings a consistent speedup @@ -283,7 +149,7 @@ def gradient(f, *varargs, **kwargs): edge_order = kwargs.pop('edge_order', 1) if kwargs: raise TypeError('"{}" are not valid keyword arguments.'.format( - '", "'.join(kwargs.keys()))) + '", "'.join(kwargs.keys()))) if edge_order > 2: raise ValueError("'edge_order' greater than 2 not supported") @@ -293,10 +159,10 @@ def gradient(f, *varargs, **kwargs): outvals = [] # create slice objects --- initially all are [:, :, ..., :] - slice1 = [slice(None)]*N - slice2 = [slice(None)]*N - slice3 = [slice(None)]*N - slice4 = [slice(None)]*N + slice1 = [slice(None)] * N + slice2 = [slice(None)] * N + slice3 = [slice(None)] * N + slice4 = [slice(None)] * N otype = f.dtype.char if otype not in ['f', 'd', 'F', 'D', 'm', 'M']: @@ -304,7 +170,8 @@ def gradient(f, *varargs, **kwargs): # Difference of datetime64 elements results in timedelta64 if otype == 'M': - # Need to use the full dtype name because it contains unit information + # Need to use the full dtype name because it contains unit + # information otype = f.dtype.name.replace('datetime', 'timedelta') elif otype == 'm': # Needs to keep the specific units, can't be a general unit @@ -321,8 +188,9 @@ def gradient(f, *varargs, **kwargs): for i, axis in enumerate(axes): if y.shape[axis] < edge_order + 1: raise ValueError( - "Shape of array too small to calculate a numerical gradient, " - "at least (edge_order + 1) elements are required.") + "Shape of array too small to calculate a numerical " + "gradient, at least (edge_order + 1) elements are " + "required.") # result allocation out = np.empty_like(y, dtype=otype) @@ -339,14 +207,15 @@ def gradient(f, *varargs, **kwargs): else: dx1 = dx[i][0:-1] dx2 = dx[i][1:] - a = -(dx2)/(dx1 * (dx1 + dx2)) + a = -(dx2) / (dx1 * (dx1 + dx2)) b = (dx2 - dx1) / (dx1 * dx2) c = dx1 / (dx2 * (dx1 + dx2)) # fix the shape for broadcasting shape = np.ones(N, dtype=int) shape[axis] = -1 a.shape = b.shape = c.shape = shape - # 1D equivalent -- out[1:-1] = a * f[:-2] + b * f[1:-1] + c * f[2:] + # 1D equivalent -- + # out[1:-1] = a * f[:-2] + b * f[1:-1] + c * f[2:] out[slice1] = a * f[slice2] + b * f[slice3] + c * f[slice4] # Numerical differentiation: 1st order edges @@ -378,7 +247,7 @@ def gradient(f, *varargs, **kwargs): else: dx1 = dx[i][0] dx2 = dx[i][1] - a = -(2. * dx1 + dx2)/(dx1 * (dx1 + dx2)) + a = -(2. * dx1 + dx2) / (dx1 * (dx1 + dx2)) b = (dx1 + dx2) / (dx1 * dx2) c = - dx1 / (dx2 * (dx1 + dx2)) # 1D equivalent -- out[0] = a * y[0] + b * y[1] + c * y[2] From d6be041aea9a72ea40659365e5565de766e3a4b8 Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Wed, 5 Sep 2018 07:29:15 +0900 Subject: [PATCH 09/24] simplify dask.gradient --- xarray/core/dask_array_ops.py | 35 ++++++++++++----------------------- 1 file changed, 12 insertions(+), 23 deletions(-) diff --git a/xarray/core/dask_array_ops.py b/xarray/core/dask_array_ops.py index 66eff43815e..bbaf0c2db9e 100644 --- a/xarray/core/dask_array_ops.py +++ b/xarray/core/dask_array_ops.py @@ -118,35 +118,24 @@ def gradient(a, coord, axis, edge_order): depth = {d: 1 if d == axis else 0 for d in range(a.ndim)} # temporary pad zero at the boundary - boundary = {d: 0 for d in range(a.ndim)} + boundary = "none" ag = overlap(a, depth=depth, boundary=boundary) n_chunk = len(a.chunks[axis]) - array_loc_stop = np.cumsum(np.array(a.chunks[axis])) + 2 + # adjust the coordinate position with overlap + array_loc_stop = np.cumsum(np.array(a.chunks[axis])) + 1 array_loc_start = array_loc_stop - np.array(a.chunks[axis]) - 2 - # extrapolate the edge point temporary - coord_pad = np.concatenate([2 * coord[: 1] - coord[1: 2], coord, - 2 * coord[-1:] - coord[-2: -1]]) - + array_loc_stop[-1] -= 1 + array_loc_start[0] = 0 + def func(x, block_id): - x = np.asarray(x) block_loc = block_id[axis] - c = coord_pad[array_loc_start[block_loc]: array_loc_stop[block_loc]] + c = coord[array_loc_start[block_loc]: array_loc_stop[block_loc]] grad = np.gradient(x, c, axis=axis, edge_order=edge_order) - - # overwrite the boundary value - if block_loc == 0: - idx = (slice(None), ) * axis + (slice(1, edge_order + 2), ) - c1 = c[1: edge_order + 2] - g1 = np.gradient(x[idx], c1, axis=axis, edge_order=edge_order) - grad[(slice(None), ) * axis + (1, )] = g1[ - (slice(None), ) * axis + (0, )] - if block_loc == n_chunk - 1: - idx = (slice(None), ) * axis + (slice(-edge_order - 2, -1), ) - c1 = c[-edge_order - 2: -1] - g1 = np.gradient(x[idx], c1, axis=axis, edge_order=edge_order) - grad[(slice(None), ) * axis + (-2, )] = g1[ - (slice(None), ) * axis + (-1, )] return grad - return trim_internal(ag.map_blocks(func, dtype=a.dtype), depth) + return a.map_overlap( + func, + dtype=a.dtype, + depth={j: 1 if j == axis else 0 for j in range(a.ndim)}, + boundary="none") From a0834600659e3223e5491cff6085838e9a86eaa4 Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Wed, 5 Sep 2018 07:36:03 +0900 Subject: [PATCH 10/24] lint --- xarray/core/dask_array_ops.py | 15 ++++----------- 1 file changed, 4 insertions(+), 11 deletions(-) diff --git a/xarray/core/dask_array_ops.py b/xarray/core/dask_array_ops.py index bbaf0c2db9e..e3afc7cb69b 100644 --- a/xarray/core/dask_array_ops.py +++ b/xarray/core/dask_array_ops.py @@ -116,18 +116,12 @@ def gradient(a, coord, axis, edge_order): raise ValueError('Chunk size must be larger than edge_order + 1. ' 'Given {}. Rechunk to proceed.'.format(c)) - depth = {d: 1 if d == axis else 0 for d in range(a.ndim)} - # temporary pad zero at the boundary - boundary = "none" - ag = overlap(a, depth=depth, boundary=boundary) - - n_chunk = len(a.chunks[axis]) # adjust the coordinate position with overlap array_loc_stop = np.cumsum(np.array(a.chunks[axis])) + 1 array_loc_start = array_loc_stop - np.array(a.chunks[axis]) - 2 array_loc_stop[-1] -= 1 array_loc_start[0] = 0 - + def func(x, block_id): block_loc = block_id[axis] c = coord[array_loc_start[block_loc]: array_loc_stop[block_loc]] @@ -135,7 +129,6 @@ def func(x, block_id): return grad return a.map_overlap( - func, - dtype=a.dtype, - depth={j: 1 if j == axis else 0 for j in range(a.ndim)}, - boundary="none") + func, dtype=a.dtype, + depth={j: 1 if j == axis else 0 for j in range(a.ndim)}, + boundary="none") From 267694db128562ef4040975e09c0257f09e17a70 Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Wed, 5 Sep 2018 10:45:12 +0900 Subject: [PATCH 11/24] Use npcompat.gradient in tests --- xarray/core/dask_array_ops.py | 4 ++-- xarray/tests/test_computation.py | 5 +++-- 2 files changed, 5 insertions(+), 4 deletions(-) diff --git a/xarray/core/dask_array_ops.py b/xarray/core/dask_array_ops.py index e3afc7cb69b..a8f2637e55b 100644 --- a/xarray/core/dask_array_ops.py +++ b/xarray/core/dask_array_ops.py @@ -3,7 +3,7 @@ import numpy as np -from . import nputils +from . import nputils, npcompat from . import dtypes try: @@ -125,7 +125,7 @@ def gradient(a, coord, axis, edge_order): def func(x, block_id): block_loc = block_id[axis] c = coord[array_loc_start[block_loc]: array_loc_stop[block_loc]] - grad = np.gradient(x, c, axis=axis, edge_order=edge_order) + grad = npcompat.gradient(x, c, axis=axis, edge_order=edge_order) return grad return a.map_overlap( diff --git a/xarray/tests/test_computation.py b/xarray/tests/test_computation.py index 364cdd0be85..e62192f562e 100644 --- a/xarray/tests/test_computation.py +++ b/xarray/tests/test_computation.py @@ -14,6 +14,7 @@ _UFuncSignature, apply_ufunc, broadcast_compat_data, collect_dict_values, join_dict_keys, ordered_set_intersection, ordered_set_union, result_name, unified_dim_sizes) +from xarray.core import npcompat from xarray.testing import assert_equal from . import raises_regex, requires_dask, has_dask @@ -994,7 +995,7 @@ def test_gradient(dask, edge_order): # along x actual = xr.gradient(da, 'x', edge_order) expected_x = xr.DataArray( - np.gradient(da, da['x'], axis=0, edge_order=edge_order), + npcompat.gradient(da, da['x'], axis=0, edge_order=edge_order), dims=da.dims, coords=da.coords) assert_equal(expected_x, actual) assert_equal(actual, ds.gradient('x', edge_order=edge_order)['var']) @@ -1004,7 +1005,7 @@ def test_gradient(dask, edge_order): # along y actual = xr.gradient(da, 'y', edge_order) expected_y = xr.DataArray( - np.gradient(da, da['y'], axis=1, edge_order=edge_order), + npcompat.gradient(da, da['y'], axis=1, edge_order=edge_order), dims=da.dims, coords=da.coords) assert_equal(expected_y, actual) assert_equal(actual, ds.gradient('y', edge_order=edge_order)['var']) From 2a71b62c4cbbe2e891febb60cb5b81f0368f477f Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Thu, 6 Sep 2018 09:24:04 +0900 Subject: [PATCH 12/24] move gradient to dask_array_compat --- xarray/core/dask_array_compat.py | 120 +++++++++++++++++++++++++++++++ xarray/core/dask_array_ops.py | 28 -------- xarray/core/duck_array_ops.py | 2 +- 3 files changed, 121 insertions(+), 29 deletions(-) diff --git a/xarray/core/dask_array_compat.py b/xarray/core/dask_array_compat.py index c2417345f55..7dc56daae0e 100644 --- a/xarray/core/dask_array_compat.py +++ b/xarray/core/dask_array_compat.py @@ -1,6 +1,9 @@ from __future__ import absolute_import, division, print_function +from distutils.version import LooseVersion + import numpy as np +from dask import __version__ as dask_version import dask.array as da try: @@ -30,3 +33,120 @@ def isin(element, test_elements, assume_unique=False, invert=False): if invert: result = ~result return result + + +if LooseVersion(dask_version) > LooseVersion('1.20'): + gradient = da.gradient + +else: # pragma: no cover + # Copied from dask v0.19.2 + # Used under the terms of Dask's license, see licenses/DASK_LICENSE. + import math + from numbers import Integral, Real + + + def validate_axis(axis, ndim): + """ Validate an input to axis= keywords """ + if isinstance(axis, (tuple, list)): + return tuple(validate_axis(ax, ndim) for ax in axis) + if not isinstance(axis, Integral): + raise TypeError("Axis value must be an integer, got %s" % axis) + if axis < -ndim or axis >= ndim: + raise AxisError("Axis %d is out of bounds for array of dimension %d" + % (axis, ndim)) + if axis < 0: + axis += ndim + return axis + + + def _gradient_kernel(x, block_id, coord, axis, array_locs, grad_kwargs): + """ + x: nd-array + array of one block + coord: 1d-array or scalar + coordinate along which the gradient is computed. + axis: int + axis along which the gradient is computed + array_locs: + actual location along axis. None if coordinate is scalar + grad_kwargs: + keyword to be passed to np.gradient + """ + block_loc = block_id[axis] + if array_locs is not None: + coord = coord[array_locs[0][block_loc]: array_locs[1][block_loc]] + grad = np.gradient(x, coord, axis=axis, **grad_kwargs) + return grad + + + def gradient(f, *varargs, **kwargs): + f = da.asarray(f) + + kwargs["edge_order"] = math.ceil(kwargs.get("edge_order", 1)) + if kwargs["edge_order"] > 2: + raise ValueError("edge_order must be less than or equal to 2.") + + drop_result_list = False + axis = kwargs.pop("axis", None) + if axis is None: + axis = tuple(range(f.ndim)) + elif isinstance(axis, Integral): + drop_result_list = True + axis = (axis,) + + axis = validate_axis(axis, f.ndim) + + if len(axis) != len(set(axis)): + raise ValueError("duplicate axes not allowed") + + axis = tuple(ax % f.ndim for ax in axis) + + if varargs == (): + varargs = (1,) + if len(varargs) == 1: + varargs = len(axis) * varargs + if len(varargs) != len(axis): + raise TypeError( + "Spacing must either be a single scalar, or a scalar / " + "1d-array per axis" + ) + + if issubclass(f.dtype.type, (np.bool8, Integral)): + f = f.astype(float) + elif issubclass(f.dtype.type, Real) and f.dtype.itemsize < 4: + f = f.astype(float) + + results = [] + for i, ax in enumerate(axis): + for c in f.chunks[ax]: + if c < kwargs["edge_order"] + 1: + raise ValueError( + 'Chunk size must be larger than edge_order + 1. ' + 'Given {}. Rechunk to proceed.'.format(c)) + + if np.isscalar(varargs[i]): + array_locs = None + else: + # coordinate position for each block taking overlap into + # account + array_loc_stop = np.cumsum(np.array(f.chunks[ax])) + 1 + array_loc_start = array_loc_stop - np.array(f.chunks[ax]) - 2 + array_loc_stop[-1] -= 1 + array_loc_start[0] = 0 + array_locs = (array_loc_start, array_loc_stop) + + results.append(f.map_overlap( + _gradient_kernel, + dtype=f.dtype, + depth={j: 1 if j == ax else 0 for j in range(f.ndim)}, + boundary="none", + coord=varargs[i], + axis=ax, + array_locs=array_locs, + grad_kwargs=kwargs, + )) + + if drop_result_list: + results = results[0] + + return results diff --git a/xarray/core/dask_array_ops.py b/xarray/core/dask_array_ops.py index a8f2637e55b..274035dc5b2 100644 --- a/xarray/core/dask_array_ops.py +++ b/xarray/core/dask_array_ops.py @@ -104,31 +104,3 @@ def func(x, window, axis=-1): index = (slice(None),) * axis + (slice(drop_size, drop_size + orig_shape[axis]), ) return out[index] - - -def gradient(a, coord, axis, edge_order): - """ Apply gradient along one dimension """ - if axis < 0: - axis = a.ndim + axis - - for c in a.chunks[axis]: - if c < edge_order + 1: - raise ValueError('Chunk size must be larger than edge_order + 1. ' - 'Given {}. Rechunk to proceed.'.format(c)) - - # adjust the coordinate position with overlap - array_loc_stop = np.cumsum(np.array(a.chunks[axis])) + 1 - array_loc_start = array_loc_stop - np.array(a.chunks[axis]) - 2 - array_loc_stop[-1] -= 1 - array_loc_start[0] = 0 - - def func(x, block_id): - block_loc = block_id[axis] - c = coord[array_loc_start[block_loc]: array_loc_stop[block_loc]] - grad = npcompat.gradient(x, c, axis=axis, edge_order=edge_order) - return grad - - return a.map_overlap( - func, dtype=a.dtype, - depth={j: 1 if j == axis else 0 for j in range(a.ndim)}, - boundary="none") diff --git a/xarray/core/duck_array_ops.py b/xarray/core/duck_array_ops.py index 0d44c4f108d..ef89dba2ab8 100644 --- a/xarray/core/duck_array_ops.py +++ b/xarray/core/duck_array_ops.py @@ -96,7 +96,7 @@ def isnull(data): def gradient(x, coord, axis, edge_order): if isinstance(x, dask_array_type): - return dask_array_ops.gradient( + return dask_array_compat.gradient( x, coord, axis=axis, edge_order=edge_order) return npcompat.gradient(x, coord, axis=axis, edge_order=edge_order) From b504da8a79bdf15e990d6b9bd4ef3e40ac24ed0f Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Thu, 6 Sep 2018 10:31:44 +0900 Subject: [PATCH 13/24] gradient -> differentiate --- doc/api.rst | 6 ++-- doc/whats-new.rst | 4 +-- xarray/__init__.py | 2 +- xarray/core/computation.py | 60 +++++++++----------------------- xarray/core/dataarray.py | 14 ++++---- xarray/core/dataset.py | 12 +++---- xarray/tests/test_computation.py | 23 +++++------- 7 files changed, 44 insertions(+), 77 deletions(-) diff --git a/doc/api.rst b/doc/api.rst index 67e396f656f..cc6dfff4e79 100644 --- a/doc/api.rst +++ b/doc/api.rst @@ -25,7 +25,7 @@ Top-level functions zeros_like ones_like dot - gradient + differentiate Dataset ======= @@ -151,7 +151,7 @@ Computation Dataset.resample Dataset.diff Dataset.quantile - Dataset.gradient + Dataset.differentiate **Aggregation**: :py:attr:`~Dataset.all` @@ -319,7 +319,7 @@ Computation DataArray.diff DataArray.dot DataArray.quantile - DataArray.gradient + DataArray.differentiate **Aggregation**: :py:attr:`~DataArray.all` diff --git a/doc/whats-new.rst b/doc/whats-new.rst index ea7176ded5e..0328f195bac 100644 --- a/doc/whats-new.rst +++ b/doc/whats-new.rst @@ -36,8 +36,8 @@ Documentation Enhancements ~~~~~~~~~~~~ -- :py:func:`~xarray.gradient`, :py:meth:`~xarray.DataArray.gradient`, and - :py:meth:`~xarray.Dataset.gradient` are newly added. +- :py:func:`~xarray.differentiate`, :py:meth:`~xarray.DataArray.differentiate`, + and :py:meth:`~xarray.Dataset.differentiate` are newly added. (:issue:`1332`) By `Keisuke Fujii `_. diff --git a/xarray/__init__.py b/xarray/__init__.py index 968b2610350..3d9d187fa80 100644 --- a/xarray/__init__.py +++ b/xarray/__init__.py @@ -10,7 +10,7 @@ from .core.alignment import align, broadcast, broadcast_arrays from .core.common import full_like, zeros_like, ones_like from .core.combine import concat, auto_combine -from .core.computation import apply_ufunc, dot, gradient, where +from .core.computation import apply_ufunc, differentiate, dot, where from .core.extensions import (register_dataarray_accessor, register_dataset_accessor) from .core.variable import as_variable, Variable, IndexVariable, Coordinate diff --git a/xarray/core/computation.py b/xarray/core/computation.py index 5a9935d145d..5744570ff35 100644 --- a/xarray/core/computation.py +++ b/xarray/core/computation.py @@ -1139,7 +1139,7 @@ def where(cond, x, y): dask='allowed') -def _gradient_once(variable, coord, edge_order): +def _differentiate_variable(variable, coord, edge_order): """ Compute the gradient along 1 dimension for Variable. variable, coord: Variable """ @@ -1151,7 +1151,7 @@ def _gradient_once(variable, coord, edge_order): return Variable(variable.dims, result_array) -def gradient(dataarray, coords, edge_order=1): +def differentiate(dataarray, coord, edge_order=1): """ Compute the gradient of the dataarray with the second order accurate central differences. @@ -1159,14 +1159,14 @@ def gradient(dataarray, coords, edge_order=1): ---------- dataarray: xr.DataArray Target array - coords: str or sequence of strs - The coordinate to be used to compute the gradient. + coord: str + The coordinate to be used to differentiate the array. edge_order: 1 or 2. Default 1 N-th order accurate differences at the boundaries. Returns ------- - gradient: DataArray or a sequence of DataArrays + differentiate: DataArray or a sequence of DataArrays See also -------- @@ -1187,7 +1187,7 @@ def gradient(dataarray, coords, edge_order=1): * x (x) float64 0.0 0.1 1.1 1.2 Dimensions without coordinates: y >>> - >>> xr.gradient(da, 'x') + >>> xr.differentiate(da, 'x') array([[30. , 30. , 30. ], [27.545455, 27.545455, 27.545455], @@ -1196,23 +1196,6 @@ def gradient(dataarray, coords, edge_order=1): Coordinates: * x (x) float64 0.0 0.1 1.1 1.2 Dimensions without coordinates: y - >>> - >>> xr.gradient(da, ('x', 'y')) - ( - array([[30. , 30. , 30. ], - [27.545455, 27.545455, 27.545455], - [27.545455, 27.545455, 27.545455], - [30. , 30. , 30. ]]) - Coordinates: - * x (x) float64 0.0 0.1 1.1 1.2 - Dimensions without coordinates: y, - array([[1., 1., 1.], - [1., 1., 1.], - [1., 1., 1.], - [1., 1., 1.]]) - Coordinates: - * x (x) float64 0.0 0.1 1.1 1.2 - Dimensions without coordinates: y) """ from .dataarray import DataArray @@ -1220,26 +1203,15 @@ def gradient(dataarray, coords, edge_order=1): raise TypeError( 'Only DataArray is supported. Given {}.'.format(type(dataarray))) - return_sequence = True - if not isinstance(coords, (tuple, list)): - coords = (coords, ) - return_sequence = False - - result = [] - for coord in coords: - if coord not in dataarray.coords and coord not in dataarray.dims: - raise ValueError('Coordiante {} does not exist.'.format(coord)) + if coord not in dataarray.coords and coord not in dataarray.dims: + raise ValueError('Coordiante {} does not exist.'.format(coord)) - coord_var = dataarray[coord].variable - if coord_var.ndim != 1: - raise ValueError( - 'Only 1d-coordinate is supported. {} is {} ' - 'dimensional.'.format(coord, dataarray[coord].ndim)) - - result.append(DataArray( - _gradient_once(dataarray.variable, coord_var, edge_order), - coords=dataarray.coords)) + coord_var = dataarray[coord].variable + if coord_var.ndim != 1: + raise ValueError( + 'Only 1d-coordinate is supported. {} is {} ' + 'dimensional.'.format(coord, dataarray[coord].ndim)) - if return_sequence: - return tuple(result) - return result[0] + return DataArray( + _differentiate_variable(dataarray.variable, coord_var, edge_order), + coords=dataarray.coords) diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py index 3defb37d2cf..5c9f78ff20b 100644 --- a/xarray/core/dataarray.py +++ b/xarray/core/dataarray.py @@ -2267,25 +2267,25 @@ def rank(self, dim, pct=False, keep_attrs=False): ds = self._to_temp_dataset().rank(dim, pct=pct, keep_attrs=keep_attrs) return self._from_temp_dataset(ds) - def gradient(self, coord, edge_order=1): - """ Compute the gradient with the second order accurate central + def differentiate(self, coord, edge_order=1): + """ Differentiate the array with the second order accurate central differences. Parameters ---------- - coords: str + coord: str The coordinate to be used to compute the gradient. edge_order: 1 or 2. Default 1 N-th order accurate differences at the boundaries. Returns ------- - gradient: DataArray + differentiated: DataArray See also -------- numpy.gradient: corresponding numpy function - xr.gradient: more numpy-like function for xarray object. + xr.differentiate: more numpy-like function for xarray object. Examples -------- @@ -2302,7 +2302,7 @@ def gradient(self, coord, edge_order=1): * x (x) float64 0.0 0.1 1.1 1.2 Dimensions without coordinates: y >>> - >>> da.gradient('x') + >>> da.differentiate('x') array([[30. , 30. , 30. ], [27.545455, 27.545455, 27.545455], @@ -2312,7 +2312,7 @@ def gradient(self, coord, edge_order=1): * x (x) float64 0.0 0.1 1.1 1.2 Dimensions without coordinates: y """ - ds = self._to_temp_dataset().gradient(coord, edge_order) + ds = self._to_temp_dataset().differentiate(coord, edge_order) return self._from_temp_dataset(ds) diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py index d7e4589360a..ae218bd05ac 100644 --- a/xarray/core/dataset.py +++ b/xarray/core/dataset.py @@ -3645,25 +3645,25 @@ def rank(self, dim, pct=False, keep_attrs=False): attrs = self.attrs if keep_attrs else None return self._replace_vars_and_dims(variables, coord_names, attrs=attrs) - def gradient(self, coord, edge_order=1): - """ Compute the gradient with the second order accurate central + def differentiate(self, coord, edge_order=1): + """ Differentiate with the second order accurate central differences. Parameters ---------- - coords: str + coord: str The coordinate to be used to compute the gradient. edge_order: 1 or 2. Default 1 N-th order accurate differences at the boundaries. Returns ------- - gradient: Dataset + differentiated: Dataset See also -------- numpy.gradient: corresponding numpy function - xr.gradient: more numpy-like function for xarray object. + xr.differentiated: more numpy-like function for xarray object. """ if coord not in self.variables and coord not in self.dims: raise ValueError('Coordinate {} does not exist.'.format(coord)) @@ -3677,7 +3677,7 @@ def gradient(self, coord, edge_order=1): variables = OrderedDict() for k, v in self.variables.items(): if k in self.data_vars and dim in v.dims: - variables[k] = computation._gradient_once( + variables[k] = computation._differentiate_variable( v, coord_var, edge_order) else: variables[k] = v diff --git a/xarray/tests/test_computation.py b/xarray/tests/test_computation.py index e62192f562e..90cc46336d1 100644 --- a/xarray/tests/test_computation.py +++ b/xarray/tests/test_computation.py @@ -993,29 +993,24 @@ def test_gradient(dask, edge_order): ds = xr.Dataset({'var': da}) # along x - actual = xr.gradient(da, 'x', edge_order) + actual = xr.differentiate(da, 'x', edge_order) expected_x = xr.DataArray( npcompat.gradient(da, da['x'], axis=0, edge_order=edge_order), dims=da.dims, coords=da.coords) assert_equal(expected_x, actual) - assert_equal(actual, ds.gradient('x', edge_order=edge_order)['var']) - assert_equal(ds['var'].gradient('x', edge_order=edge_order), - ds.gradient('x', edge_order=edge_order)['var']) + assert_equal(actual, ds.differentiate('x', edge_order=edge_order)['var']) + assert_equal(ds['var'].differentiate('x', edge_order=edge_order), + ds.differentiate('x', edge_order=edge_order)['var']) # along y - actual = xr.gradient(da, 'y', edge_order) + actual = xr.differentiate(da, 'y', edge_order) expected_y = xr.DataArray( npcompat.gradient(da, da['y'], axis=1, edge_order=edge_order), dims=da.dims, coords=da.coords) assert_equal(expected_y, actual) - assert_equal(actual, ds.gradient('y', edge_order=edge_order)['var']) - assert_equal(ds['var'].gradient('y', edge_order=edge_order), - ds.gradient('y', edge_order=edge_order)['var']) - - # along x and y - actual = xr.gradient(da, ('x', 'y'), edge_order) - assert_equal(expected_x, actual[0]) - assert_equal(expected_y, actual[1]) + assert_equal(actual, ds.differentiate('y', edge_order=edge_order)['var']) + assert_equal(ds['var'].differentiate('y', edge_order=edge_order), + ds.differentiate('y', edge_order=edge_order)['var']) with pytest.raises(ValueError): - xr.gradient(da, ('x2d')) + xr.differentiate(da, ('x2d')) From fb356c5e58bdb3851eef5148056f261bfd146d37 Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Thu, 6 Sep 2018 11:06:47 +0900 Subject: [PATCH 14/24] lint --- xarray/core/dask_array_compat.py | 7 +++---- xarray/core/dask_array_ops.py | 2 +- 2 files changed, 4 insertions(+), 5 deletions(-) diff --git a/xarray/core/dask_array_compat.py b/xarray/core/dask_array_compat.py index 7dc56daae0e..c6e5717adfb 100644 --- a/xarray/core/dask_array_compat.py +++ b/xarray/core/dask_array_compat.py @@ -44,6 +44,7 @@ def isin(element, test_elements, assume_unique=False, invert=False): import math from numbers import Integral, Real + AxisError = np.AxisError def validate_axis(axis, ndim): """ Validate an input to axis= keywords """ @@ -52,13 +53,12 @@ def validate_axis(axis, ndim): if not isinstance(axis, Integral): raise TypeError("Axis value must be an integer, got %s" % axis) if axis < -ndim or axis >= ndim: - raise AxisError("Axis %d is out of bounds for array of dimension %d" - % (axis, ndim)) + raise AxisError("Axis %d is out of bounds for array of dimension " + "%d" % (axis, ndim)) if axis < 0: axis += ndim return axis - def _gradient_kernel(x, block_id, coord, axis, array_locs, grad_kwargs): """ x: nd-array @@ -78,7 +78,6 @@ def _gradient_kernel(x, block_id, coord, axis, array_locs, grad_kwargs): grad = np.gradient(x, coord, axis=axis, **grad_kwargs) return grad - def gradient(f, *varargs, **kwargs): f = da.asarray(f) diff --git a/xarray/core/dask_array_ops.py b/xarray/core/dask_array_ops.py index 274035dc5b2..423a65aa3c2 100644 --- a/xarray/core/dask_array_ops.py +++ b/xarray/core/dask_array_ops.py @@ -3,7 +3,7 @@ import numpy as np -from . import nputils, npcompat +from . import nputils from . import dtypes try: From 7a0b57ff24bf01d7a2d376092cc5dca05d4995be Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Mon, 10 Sep 2018 14:44:56 +0900 Subject: [PATCH 15/24] Update dask_array_compat --- xarray/core/dask_array_compat.py | 17 +++++++++++------ 1 file changed, 11 insertions(+), 6 deletions(-) diff --git a/xarray/core/dask_array_compat.py b/xarray/core/dask_array_compat.py index c6e5717adfb..5e6b81a253d 100644 --- a/xarray/core/dask_array_compat.py +++ b/xarray/core/dask_array_compat.py @@ -35,7 +35,7 @@ def isin(element, test_elements, assume_unique=False, invert=False): return result -if LooseVersion(dask_version) > LooseVersion('1.20'): +if LooseVersion(dask_version) > LooseVersion('1.19.2'): gradient = da.gradient else: # pragma: no cover @@ -74,7 +74,7 @@ def _gradient_kernel(x, block_id, coord, axis, array_locs, grad_kwargs): """ block_loc = block_id[axis] if array_locs is not None: - coord = coord[array_locs[0][block_loc]: array_locs[1][block_loc]] + coord = coord[array_locs[0][block_loc]:array_locs[1][block_loc]] grad = np.gradient(x, coord, axis=axis, **grad_kwargs) return grad @@ -118,18 +118,23 @@ def gradient(f, *varargs, **kwargs): results = [] for i, ax in enumerate(axis): for c in f.chunks[ax]: - if c < kwargs["edge_order"] + 1: + if np.min(c) < kwargs["edge_order"] + 1: raise ValueError( 'Chunk size must be larger than edge_order + 1. ' - 'Given {}. Rechunk to proceed.'.format(c)) + 'Minimum chunk for aixs {} is {}. Rechunk to ' + 'proceed.'.format(np.min(c), ax)) if np.isscalar(varargs[i]): array_locs = None else: + if isinstance(varargs[i], da.Array): + raise NotImplementedError( + 'dask array coordinated is not supported.') # coordinate position for each block taking overlap into # account - array_loc_stop = np.cumsum(np.array(f.chunks[ax])) + 1 - array_loc_start = array_loc_stop - np.array(f.chunks[ax]) - 2 + chunk = np.array(f.chunks[ax]) + array_loc_stop = np.cumsum(chunk) + 1 + array_loc_start = array_loc_stop - chunk - 2 array_loc_stop[-1] -= 1 array_loc_start[0] = 0 array_locs = (array_loc_start, array_loc_stop) From e93b926c14bdb6a17ab656bffad48eb659e6c422 Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Mon, 10 Sep 2018 15:02:55 +0900 Subject: [PATCH 16/24] Added a link from diff --- xarray/core/dataarray.py | 3 +++ xarray/core/dataset.py | 3 +++ 2 files changed, 6 insertions(+) diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py index 4bf1f09a796..a31b85592c7 100644 --- a/xarray/core/dataarray.py +++ b/xarray/core/dataarray.py @@ -2010,6 +2010,9 @@ def diff(self, dim, n=1, label='upper'): Coordinates: * x (x) int64 3 4 + See Also + -------- + DataArray.differentiate """ ds = self._to_temp_dataset().diff(n=n, dim=dim, label=label) return self._from_temp_dataset(ds) diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py index f4a18cb5ef8..3058cb1aaa1 100644 --- a/xarray/core/dataset.py +++ b/xarray/core/dataset.py @@ -3313,6 +3313,9 @@ def diff(self, dim, n=1, label='upper'): Data variables: foo (x) int64 1 -1 + See Also + -------- + Dataset.differentiate """ if n == 0: return self From 74acc81e1046e8f91f10477e95efb4c05eea3f4c Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Wed, 12 Sep 2018 11:29:59 +0900 Subject: [PATCH 17/24] remove xr.differentiate --- doc/api.rst | 1 - xarray/__init__.py | 2 +- xarray/core/computation.py | 78 -------------------------------- xarray/core/dataarray.py | 1 - xarray/core/dataset.py | 13 ++++-- xarray/tests/test_computation.py | 38 ---------------- xarray/tests/test_dataset.py | 37 ++++++++++++++- 7 files changed, 45 insertions(+), 125 deletions(-) diff --git a/doc/api.rst b/doc/api.rst index cc6dfff4e79..ce15e8f9d49 100644 --- a/doc/api.rst +++ b/doc/api.rst @@ -25,7 +25,6 @@ Top-level functions zeros_like ones_like dot - differentiate Dataset ======= diff --git a/xarray/__init__.py b/xarray/__init__.py index 3d9d187fa80..27a3f6b1eff 100644 --- a/xarray/__init__.py +++ b/xarray/__init__.py @@ -10,7 +10,7 @@ from .core.alignment import align, broadcast, broadcast_arrays from .core.common import full_like, zeros_like, ones_like from .core.combine import concat, auto_combine -from .core.computation import apply_ufunc, differentiate, dot, where +from .core.computation import apply_ufunc, dot, where from .core.extensions import (register_dataarray_accessor, register_dataset_accessor) from .core.variable import as_variable, Variable, IndexVariable, Coordinate diff --git a/xarray/core/computation.py b/xarray/core/computation.py index 5744570ff35..bdba72cb48a 100644 --- a/xarray/core/computation.py +++ b/xarray/core/computation.py @@ -1137,81 +1137,3 @@ def where(cond, x, y): join='exact', dataset_join='exact', dask='allowed') - - -def _differentiate_variable(variable, coord, edge_order): - """ Compute the gradient along 1 dimension for Variable. - variable, coord: Variable - """ - from .variable import Variable - - result_array = duck_array_ops.gradient( - variable.data, coord.data, edge_order=edge_order, - axis=variable.get_axis_num(coord.dims[0])) - return Variable(variable.dims, result_array) - - -def differentiate(dataarray, coord, edge_order=1): - """ Compute the gradient of the dataarray with the second order accurate - central differences. - - Parameters - ---------- - dataarray: xr.DataArray - Target array - coord: str - The coordinate to be used to differentiate the array. - edge_order: 1 or 2. Default 1 - N-th order accurate differences at the boundaries. - - Returns - ------- - differentiate: DataArray or a sequence of DataArrays - - See also - -------- - numpy.gradient: corresponding numpy function - - Examples - -------- - - >>> da = xr.DataArray(np.arange(12).reshape(4, 3), dims=['x', 'y'], - ... coords={'x': [0, 0.1, 1.1, 1.2]}) - >>> da - - array([[ 0, 1, 2], - [ 3, 4, 5], - [ 6, 7, 8], - [ 9, 10, 11]]) - Coordinates: - * x (x) float64 0.0 0.1 1.1 1.2 - Dimensions without coordinates: y - >>> - >>> xr.differentiate(da, 'x') - - array([[30. , 30. , 30. ], - [27.545455, 27.545455, 27.545455], - [27.545455, 27.545455, 27.545455], - [30. , 30. , 30. ]]) - Coordinates: - * x (x) float64 0.0 0.1 1.1 1.2 - Dimensions without coordinates: y - """ - from .dataarray import DataArray - - if not isinstance(dataarray, DataArray): - raise TypeError( - 'Only DataArray is supported. Given {}.'.format(type(dataarray))) - - if coord not in dataarray.coords and coord not in dataarray.dims: - raise ValueError('Coordiante {} does not exist.'.format(coord)) - - coord_var = dataarray[coord].variable - if coord_var.ndim != 1: - raise ValueError( - 'Only 1d-coordinate is supported. {} is {} ' - 'dimensional.'.format(coord, dataarray[coord].ndim)) - - return DataArray( - _differentiate_variable(dataarray.variable, coord_var, edge_order), - coords=dataarray.coords) diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py index a31b85592c7..33be8a5f4ab 100644 --- a/xarray/core/dataarray.py +++ b/xarray/core/dataarray.py @@ -2310,7 +2310,6 @@ def differentiate(self, coord, edge_order=1): See also -------- numpy.gradient: corresponding numpy function - xr.differentiate: more numpy-like function for xarray object. Examples -------- diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py index 3058cb1aaa1..3070282670e 100644 --- a/xarray/core/dataset.py +++ b/xarray/core/dataset.py @@ -3684,22 +3684,25 @@ def differentiate(self, coord, edge_order=1): See also -------- numpy.gradient: corresponding numpy function - xr.differentiated: more numpy-like function for xarray object. """ + from .variable import Variable + if coord not in self.variables and coord not in self.dims: raise ValueError('Coordinate {} does not exist.'.format(coord)) coord_var = self[coord].variable if coord_var.ndim != 1: - raise ValueError('Coordinate {} must be 1 dimensional but {} is {}' - ' dimensional'.format(coord, coord_var.ndims)) + raise ValueError('Coordinate {} must be 1 dimensional but is {}' + ' dimensional'.format(coord, coord_var.ndim)) dim = coord_var.dims[0] variables = OrderedDict() for k, v in self.variables.items(): if k in self.data_vars and dim in v.dims: - variables[k] = computation._differentiate_variable( - v, coord_var, edge_order) + grad = duck_array_ops.gradient( + v.data, coord_var.data, edge_order=edge_order, + axis=v.get_axis_num(dim)) + variables[k] = Variable(v.dims, grad) else: variables[k] = v return self._replace_vars_and_dims(variables) diff --git a/xarray/tests/test_computation.py b/xarray/tests/test_computation.py index 90cc46336d1..ca8e4e59737 100644 --- a/xarray/tests/test_computation.py +++ b/xarray/tests/test_computation.py @@ -14,8 +14,6 @@ _UFuncSignature, apply_ufunc, broadcast_compat_data, collect_dict_values, join_dict_keys, ordered_set_intersection, ordered_set_union, result_name, unified_dim_sizes) -from xarray.core import npcompat -from xarray.testing import assert_equal from . import raises_regex, requires_dask, has_dask @@ -978,39 +976,3 @@ def test_where(): actual = xr.where(cond, 1, 0) expected = xr.DataArray([1, 0], dims='x') assert_identical(expected, actual) - - -@pytest.mark.parametrize('dask', [True, False]) -@pytest.mark.parametrize('edge_order', [1, 2]) -def test_gradient(dask, edge_order): - rs = np.random.RandomState(42) - da = xr.DataArray(rs.randn(8, 6), dims=['x', 'y'], - coords={'x': [0.2, 0.35, 0.4, 0.6, 0.7, 0.75, 0.76, 0.8], - 'z': 3, 'x2d': (('x', 'y'), rs.randn(8, 6))}) - if dask and has_dask: - da = da.chunk({'x': 4}) - - ds = xr.Dataset({'var': da}) - - # along x - actual = xr.differentiate(da, 'x', edge_order) - expected_x = xr.DataArray( - npcompat.gradient(da, da['x'], axis=0, edge_order=edge_order), - dims=da.dims, coords=da.coords) - assert_equal(expected_x, actual) - assert_equal(actual, ds.differentiate('x', edge_order=edge_order)['var']) - assert_equal(ds['var'].differentiate('x', edge_order=edge_order), - ds.differentiate('x', edge_order=edge_order)['var']) - - # along y - actual = xr.differentiate(da, 'y', edge_order) - expected_y = xr.DataArray( - npcompat.gradient(da, da['y'], axis=1, edge_order=edge_order), - dims=da.dims, coords=da.coords) - assert_equal(expected_y, actual) - assert_equal(actual, ds.differentiate('y', edge_order=edge_order)['var']) - assert_equal(ds['var'].differentiate('y', edge_order=edge_order), - ds.differentiate('y', edge_order=edge_order)['var']) - - with pytest.raises(ValueError): - xr.differentiate(da, ('x2d')) diff --git a/xarray/tests/test_dataset.py b/xarray/tests/test_dataset.py index d22d8470dc6..bd8a0ab1666 100644 --- a/xarray/tests/test_dataset.py +++ b/xarray/tests/test_dataset.py @@ -15,7 +15,7 @@ from xarray import ( DataArray, Dataset, IndexVariable, MergeError, Variable, align, backends, broadcast, open_dataset, set_options) -from xarray.core import indexing, utils +from xarray.core import indexing, npcompat, utils from xarray.core.common import full_like from xarray.core.pycompat import ( OrderedDict, integer_types, iteritems, unicode_type) @@ -4492,3 +4492,38 @@ def test_raise_no_warning_for_nan_in_binary_ops(): with pytest.warns(None) as record: Dataset(data_vars={'x': ('y', [1, 2, np.NaN])}) > 0 assert len(record) == 0 + + +@pytest.mark.parametrize('dask', [True, False]) +@pytest.mark.parametrize('edge_order', [1, 2]) +def test_gradient(dask, edge_order): + rs = np.random.RandomState(42) + da = xr.DataArray(rs.randn(8, 6), dims=['x', 'y'], + coords={'x': [0.2, 0.35, 0.4, 0.6, 0.7, 0.75, 0.76, 0.8], + 'z': 3, 'x2d': (('x', 'y'), rs.randn(8, 6))}) + if dask and has_dask: + da = da.chunk({'x': 4}) + + ds = xr.Dataset({'var': da}) + + # along x + actual = da.differentiate('x', edge_order) + expected_x = xr.DataArray( + npcompat.gradient(da, da['x'], axis=0, edge_order=edge_order), + dims=da.dims, coords=da.coords) + assert_equal(expected_x, actual) + assert_equal(ds['var'].differentiate('x', edge_order=edge_order), + ds.differentiate('x', edge_order=edge_order)['var']) + + # along y + actual = da.differentiate('y', edge_order) + expected_y = xr.DataArray( + npcompat.gradient(da, da['y'], axis=1, edge_order=edge_order), + dims=da.dims, coords=da.coords) + assert_equal(expected_y, actual) + assert_equal(actual, ds.differentiate('y', edge_order=edge_order)['var']) + assert_equal(ds['var'].differentiate('y', edge_order=edge_order), + ds.differentiate('y', edge_order=edge_order)['var']) + + with pytest.raises(ValueError): + da.differentiate('x2d') From 8817306801d99538502bec78ba816cb7698ced2e Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Wed, 12 Sep 2018 23:30:42 +0900 Subject: [PATCH 18/24] Added datetime support --- xarray/core/dataarray.py | 5 +++-- xarray/core/dataset.py | 21 +++++++++++++++---- xarray/core/utils.py | 21 +++++++++++++++++++ xarray/tests/test_dataset.py | 40 +++++++++++++++++++++++++++++++++++- 4 files changed, 80 insertions(+), 7 deletions(-) diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py index 33be8a5f4ab..0fc25b5a051 100644 --- a/xarray/core/dataarray.py +++ b/xarray/core/dataarray.py @@ -2292,7 +2292,7 @@ def rank(self, dim, pct=False, keep_attrs=False): ds = self._to_temp_dataset().rank(dim, pct=pct, keep_attrs=keep_attrs) return self._from_temp_dataset(ds) - def differentiate(self, coord, edge_order=1): + def differentiate(self, coord, edge_order=1, time_unit=None): """ Differentiate the array with the second order accurate central differences. @@ -2336,7 +2336,8 @@ def differentiate(self, coord, edge_order=1): * x (x) float64 0.0 0.1 1.1 1.2 Dimensions without coordinates: y """ - ds = self._to_temp_dataset().differentiate(coord, edge_order) + ds = self._to_temp_dataset().differentiate( + coord, edge_order, time_unit) return self._from_temp_dataset(ds) diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py index 3070282670e..0fa751ed58f 100644 --- a/xarray/core/dataset.py +++ b/xarray/core/dataset.py @@ -31,7 +31,7 @@ OrderedDict, basestring, dask_array_type, integer_types, iteritems, range) from .utils import ( Frozen, SortedKeysDict, either_dict_or_kwargs, decode_numpy_dict_values, - ensure_us_time_resolution, hashable, maybe_wrap_array) + ensure_us_time_resolution, hashable, maybe_wrap_array, to_numeric) from .variable import IndexVariable, Variable, as_variable, broadcast_variables # list of attributes of pd.DatetimeIndex that are ndarrays of time info @@ -3666,7 +3666,7 @@ def rank(self, dim, pct=False, keep_attrs=False): attrs = self.attrs if keep_attrs else None return self._replace_vars_and_dims(variables, coord_names, attrs=attrs) - def differentiate(self, coord, edge_order=1): + def differentiate(self, coord, edge_order=1, time_unit=None): """ Differentiate with the second order accurate central differences. @@ -3676,6 +3676,9 @@ def differentiate(self, coord, edge_order=1): The coordinate to be used to compute the gradient. edge_order: 1 or 2. Default 1 N-th order accurate differences at the boundaries. + time_unit: None or any of {'Y', 'M', 'W', 'D', 'h', 'm', 's', 'ms', + 'us', 'ns', 'ps', 'fs', 'as'} + Unit to compute gradient. Only valid for datetime coordinate. Returns ------- @@ -3696,11 +3699,21 @@ def differentiate(self, coord, edge_order=1): ' dimensional'.format(coord, coord_var.ndim)) dim = coord_var.dims[0] + coord_data = coord_var.data + if coord_data.dtype.kind in ['m', 'M']: + if time_unit is None: + time_unit = np.datetime_data(coord_data.dtype)[0] + coord_data = to_numeric(coord_data, offset=True, + time_unit=time_unit) + variables = OrderedDict() for k, v in self.variables.items(): - if k in self.data_vars and dim in v.dims: + if (k in self.data_vars and dim in v.dims and + k not in self.coords): + if v.dtype.kind in ['m', 'M']: + v = to_numeric(v, offset=True, time_unit=time_unit) grad = duck_array_ops.gradient( - v.data, coord_var.data, edge_order=edge_order, + v.data, coord_data, edge_order=edge_order, axis=v.get_axis_num(dim)) variables[k] = Variable(v.dims, grad) else: diff --git a/xarray/core/utils.py b/xarray/core/utils.py index c3bb747fac5..46aa67d41df 100644 --- a/xarray/core/utils.py +++ b/xarray/core/utils.py @@ -591,3 +591,24 @@ def __iter__(self): def __len__(self): num_hidden = sum([k in self._hidden_keys for k in self._data]) return len(self._data) - num_hidden + + +def to_numeric(array, offset=False, time_unit=None, dtype=float): + """ + Make datetime array float + + offset: True to subtract minimum values to reduce round off error + time_unit: None or any of {'Y', 'M', 'W', 'D', 'h', 'm', 's', 'ms', + 'us', 'ns', 'ps', 'fs', 'as'} + dtype: target dtype + """ + if array.dtype.kind not in ['m', 'M']: + return array.astype(dtype) + if offset: + array = array - np.min(array) + else: + array = array - np.zeros(array.shape, dtype=array.dtype) + + if time_unit: + return (array / np.timedelta64(1, time_unit)).astype(dtype) + return array.astype(dtype) diff --git a/xarray/tests/test_dataset.py b/xarray/tests/test_dataset.py index bd8a0ab1666..20145b2729c 100644 --- a/xarray/tests/test_dataset.py +++ b/xarray/tests/test_dataset.py @@ -4498,8 +4498,10 @@ def test_raise_no_warning_for_nan_in_binary_ops(): @pytest.mark.parametrize('edge_order', [1, 2]) def test_gradient(dask, edge_order): rs = np.random.RandomState(42) + coord = [0.2, 0.35, 0.4, 0.6, 0.7, 0.75, 0.76, 0.8] + da = xr.DataArray(rs.randn(8, 6), dims=['x', 'y'], - coords={'x': [0.2, 0.35, 0.4, 0.6, 0.7, 0.75, 0.76, 0.8], + coords={'x': coord, 'z': 3, 'x2d': (('x', 'y'), rs.randn(8, 6))}) if dask and has_dask: da = da.chunk({'x': 4}) @@ -4514,6 +4516,8 @@ def test_gradient(dask, edge_order): assert_equal(expected_x, actual) assert_equal(ds['var'].differentiate('x', edge_order=edge_order), ds.differentiate('x', edge_order=edge_order)['var']) + # coordinate should not change + assert_equal(da['x'], actual['x']) # along y actual = da.differentiate('y', edge_order) @@ -4527,3 +4531,37 @@ def test_gradient(dask, edge_order): with pytest.raises(ValueError): da.differentiate('x2d') + + +@pytest.mark.parametrize('dask', [True, False]) +def test_gradient_datetime(dask): + rs = np.random.RandomState(42) + coord = np.array( + ['2004-07-13', '2006-01-13', '2010-08-13', '2010-09-13', + '2010-10-11', '2010-12-13', '2011-02-13', '2012-08-13'], + dtype='datetime64') + + da = xr.DataArray(rs.randn(8, 6), dims=['x', 'y'], + coords={'x': coord, + 'z': 3, 'x2d': (('x', 'y'), rs.randn(8, 6))}) + if dask and has_dask: + da = da.chunk({'x': 4}) + + # along x + actual = da.differentiate('x', edge_order=1, time_unit='D') + expected_x = xr.DataArray( + npcompat.gradient( + da, utils.to_numeric( + da['x'], offset=True, time_unit='D'), axis=0, edge_order=1), + dims=da.dims, coords=da.coords) + assert_equal(expected_x, actual) + + # for datetime variable + actual = da['x'].differentiate('x', edge_order=1, time_unit='D') + assert np.allclose(actual, 1.0) + + # with different date unit + da = xr.DataArray(coord.astype('datetime64[ms]'), dims=['x'], + coords={'x': coord}) + actual = da.differentiate('x', edge_order=1) + assert np.allclose(actual, 1.0) From 4112cd94abafd67f4c5110136526239892f9fc6d Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Thu, 13 Sep 2018 07:09:55 +0900 Subject: [PATCH 19/24] Update via comment. Use utils.to_numeric also in interp --- xarray/core/dataarray.py | 3 +++ xarray/core/dataset.py | 9 ++++----- xarray/core/missing.py | 6 +++--- xarray/core/utils.py | 14 +++++++------- xarray/tests/test_dataset.py | 8 +++++--- 5 files changed, 22 insertions(+), 18 deletions(-) diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py index 0fc25b5a051..49004ae583d 100644 --- a/xarray/core/dataarray.py +++ b/xarray/core/dataarray.py @@ -2302,6 +2302,9 @@ def differentiate(self, coord, edge_order=1, time_unit=None): The coordinate to be used to compute the gradient. edge_order: 1 or 2. Default 1 N-th order accurate differences at the boundaries. + time_unit: None or any of {'Y', 'M', 'W', 'D', 'h', 'm', 's', 'ms', + 'us', 'ns', 'ps', 'fs', 'as'} + Unit to compute gradient. Only valid for datetime coordinate. Returns ------- diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py index 0fa751ed58f..4ae530c09ea 100644 --- a/xarray/core/dataset.py +++ b/xarray/core/dataset.py @@ -3702,16 +3702,15 @@ def differentiate(self, coord, edge_order=1, time_unit=None): coord_data = coord_var.data if coord_data.dtype.kind in ['m', 'M']: if time_unit is None: - time_unit = np.datetime_data(coord_data.dtype)[0] - coord_data = to_numeric(coord_data, offset=True, - time_unit=time_unit) + time_unit, _ = np.datetime_data(coord_data.dtype) + coord_data = to_numeric(coord_data, time_unit=time_unit) variables = OrderedDict() for k, v in self.variables.items(): if (k in self.data_vars and dim in v.dims and k not in self.coords): - if v.dtype.kind in ['m', 'M']: - v = to_numeric(v, offset=True, time_unit=time_unit) + if v.dtype.kind in 'mM': + v = to_numeric(v, time_unit=time_unit) grad = duck_array_ops.gradient( v.data, coord_data, edge_order=edge_order, axis=v.get_axis_num(dim)) diff --git a/xarray/core/missing.py b/xarray/core/missing.py index 90aa4ffaeda..afb34d99115 100644 --- a/xarray/core/missing.py +++ b/xarray/core/missing.py @@ -11,7 +11,7 @@ from . import rolling from .computation import apply_ufunc from .pycompat import iteritems -from .utils import is_scalar, OrderedSet +from .utils import is_scalar, OrderedSet, to_numeric from .variable import Variable, broadcast_variables from .duck_array_ops import dask_array_type @@ -414,8 +414,8 @@ def _floatize_x(x, new_x): # offset (min(x)) and the variation (x - min(x)) can be # represented by float. xmin = np.min(x[i]) - x[i] = (x[i] - xmin).astype(np.float64) - new_x[i] = (new_x[i] - xmin).astype(np.float64) + x[i] = to_numeric(x[i], offset=xmin, dtype=np.float64) + new_x[i] = to_numeric(new_x[i], offset=xmin, dtype=np.float64) return x, new_x diff --git a/xarray/core/utils.py b/xarray/core/utils.py index 46aa67d41df..dd8bfd0f6fc 100644 --- a/xarray/core/utils.py +++ b/xarray/core/utils.py @@ -593,22 +593,22 @@ def __len__(self): return len(self._data) - num_hidden -def to_numeric(array, offset=False, time_unit=None, dtype=float): +def to_numeric(array, offset=None, time_unit=None, dtype=float): """ Make datetime array float - offset: True to subtract minimum values to reduce round off error + offset: Scalar with the same type of array or None + If None, subtract minimum values to reduce round off error time_unit: None or any of {'Y', 'M', 'W', 'D', 'h', 'm', 's', 'ms', 'us', 'ns', 'ps', 'fs', 'as'} dtype: target dtype """ if array.dtype.kind not in ['m', 'M']: return array.astype(dtype) - if offset: - array = array - np.min(array) - else: - array = array - np.zeros(array.shape, dtype=array.dtype) - + if offset is None: + offset = np.min(array) + array = array - offset + if time_unit: return (array / np.timedelta64(1, time_unit)).astype(dtype) return array.astype(dtype) diff --git a/xarray/tests/test_dataset.py b/xarray/tests/test_dataset.py index 20145b2729c..f315fee92cc 100644 --- a/xarray/tests/test_dataset.py +++ b/xarray/tests/test_dataset.py @@ -4551,11 +4551,13 @@ def test_gradient_datetime(dask): actual = da.differentiate('x', edge_order=1, time_unit='D') expected_x = xr.DataArray( npcompat.gradient( - da, utils.to_numeric( - da['x'], offset=True, time_unit='D'), axis=0, edge_order=1), - dims=da.dims, coords=da.coords) + da, utils.to_numeric(da['x'], time_unit='D'), + axis=0, edge_order=1), dims=da.dims, coords=da.coords) assert_equal(expected_x, actual) + actual2 = da.differentiate('x', edge_order=1, time_unit='h') + assert np.allclose(actual, actual2 * 24) + # for datetime variable actual = da['x'].differentiate('x', edge_order=1, time_unit='D') assert np.allclose(actual, 1.0) From 1c2c88c85c1eedebe1e1ab67a9be642e38d95289 Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Tue, 18 Sep 2018 08:37:27 +0900 Subject: [PATCH 20/24] time_unit -> datetime_unit --- xarray/core/dataarray.py | 6 +++--- xarray/core/dataset.py | 14 +++++++------- xarray/core/utils.py | 8 ++++---- 3 files changed, 14 insertions(+), 14 deletions(-) diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py index 49004ae583d..8f95912704e 100644 --- a/xarray/core/dataarray.py +++ b/xarray/core/dataarray.py @@ -2292,7 +2292,7 @@ def rank(self, dim, pct=False, keep_attrs=False): ds = self._to_temp_dataset().rank(dim, pct=pct, keep_attrs=keep_attrs) return self._from_temp_dataset(ds) - def differentiate(self, coord, edge_order=1, time_unit=None): + def differentiate(self, coord, edge_order=1, datetime_unit=None): """ Differentiate the array with the second order accurate central differences. @@ -2302,7 +2302,7 @@ def differentiate(self, coord, edge_order=1, time_unit=None): The coordinate to be used to compute the gradient. edge_order: 1 or 2. Default 1 N-th order accurate differences at the boundaries. - time_unit: None or any of {'Y', 'M', 'W', 'D', 'h', 'm', 's', 'ms', + datetime_unit: None or any of {'Y', 'M', 'W', 'D', 'h', 'm', 's', 'ms', 'us', 'ns', 'ps', 'fs', 'as'} Unit to compute gradient. Only valid for datetime coordinate. @@ -2340,7 +2340,7 @@ def differentiate(self, coord, edge_order=1, time_unit=None): Dimensions without coordinates: y """ ds = self._to_temp_dataset().differentiate( - coord, edge_order, time_unit) + coord, edge_order, datetime_unit) return self._from_temp_dataset(ds) diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py index 4ae530c09ea..5b9044a4a28 100644 --- a/xarray/core/dataset.py +++ b/xarray/core/dataset.py @@ -3666,7 +3666,7 @@ def rank(self, dim, pct=False, keep_attrs=False): attrs = self.attrs if keep_attrs else None return self._replace_vars_and_dims(variables, coord_names, attrs=attrs) - def differentiate(self, coord, edge_order=1, time_unit=None): + def differentiate(self, coord, edge_order=1, datetime_unit=None): """ Differentiate with the second order accurate central differences. @@ -3676,7 +3676,7 @@ def differentiate(self, coord, edge_order=1, time_unit=None): The coordinate to be used to compute the gradient. edge_order: 1 or 2. Default 1 N-th order accurate differences at the boundaries. - time_unit: None or any of {'Y', 'M', 'W', 'D', 'h', 'm', 's', 'ms', + datetime_unit: None or any of {'Y', 'M', 'W', 'D', 'h', 'm', 's', 'ms', 'us', 'ns', 'ps', 'fs', 'as'} Unit to compute gradient. Only valid for datetime coordinate. @@ -3700,17 +3700,17 @@ def differentiate(self, coord, edge_order=1, time_unit=None): dim = coord_var.dims[0] coord_data = coord_var.data - if coord_data.dtype.kind in ['m', 'M']: - if time_unit is None: - time_unit, _ = np.datetime_data(coord_data.dtype) - coord_data = to_numeric(coord_data, time_unit=time_unit) + if coord_data.dtype.kind in 'mM': + if datetime_unit is None: + datetime_unit, _ = np.datetime_data(coord_data.dtype) + coord_data = to_numeric(coord_data, datetime_unit=datetime_unit) variables = OrderedDict() for k, v in self.variables.items(): if (k in self.data_vars and dim in v.dims and k not in self.coords): if v.dtype.kind in 'mM': - v = to_numeric(v, time_unit=time_unit) + v = to_numeric(v, datetime_unit=datetime_unit) grad = duck_array_ops.gradient( v.data, coord_data, edge_order=edge_order, axis=v.get_axis_num(dim)) diff --git a/xarray/core/utils.py b/xarray/core/utils.py index dd8bfd0f6fc..9d129d5c4f4 100644 --- a/xarray/core/utils.py +++ b/xarray/core/utils.py @@ -593,13 +593,13 @@ def __len__(self): return len(self._data) - num_hidden -def to_numeric(array, offset=None, time_unit=None, dtype=float): +def to_numeric(array, offset=None, datetime_unit=None, dtype=float): """ Make datetime array float offset: Scalar with the same type of array or None If None, subtract minimum values to reduce round off error - time_unit: None or any of {'Y', 'M', 'W', 'D', 'h', 'm', 's', 'ms', + datetime_unit: None or any of {'Y', 'M', 'W', 'D', 'h', 'm', 's', 'ms', 'us', 'ns', 'ps', 'fs', 'as'} dtype: target dtype """ @@ -609,6 +609,6 @@ def to_numeric(array, offset=None, time_unit=None, dtype=float): offset = np.min(array) array = array - offset - if time_unit: - return (array / np.timedelta64(1, time_unit)).astype(dtype) + if datetime_unit: + return (array / np.timedelta64(1, datetime_unit)).astype(dtype) return array.astype(dtype) From cbfecb4df922519aec5da9c4ffc3a6ba42dff371 Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Tue, 18 Sep 2018 08:47:32 +0900 Subject: [PATCH 21/24] Some more info in docs. --- doc/computation.rst | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/doc/computation.rst b/doc/computation.rst index 6793e667e06..125c585edb0 100644 --- a/doc/computation.rst +++ b/doc/computation.rst @@ -200,6 +200,26 @@ You can also use ``construct`` to compute a weighted rolling sum: To avoid this, use ``skipna=False`` as the above example. +Computation using Coordinates +============================= + +Xarray objects have some handy methods for the computation with their +coordinates. :py:meth:`~xarray.DataArray.differentiate` computes derivatives by +finite central differences using their coordinates, + +.. ipython:: python + a = xr.DataArray([0, 1, 2, 3], dims=['x'], coords=[0.1, 0.11, 0.2, 0.3]) + a + a.differentiate('x') + +This method can be used also for multidimensional arrays, + +.. ipython:: python + a = xr.DataArray(np.arange(8).reshape(4, 2), dims=['x'], + coords=[0.1, 0.11, 0.2, 0.3]) + a.differentiate('x') + + .. _compute.broadcasting: Broadcasting by dimension name From 2e8db19f5dedac43ad681d69d5a17441ea062a7b Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Tue, 18 Sep 2018 20:52:11 +0900 Subject: [PATCH 22/24] update test --- xarray/tests/test_dataset.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/xarray/tests/test_dataset.py b/xarray/tests/test_dataset.py index f315fee92cc..5ea6b887b5e 100644 --- a/xarray/tests/test_dataset.py +++ b/xarray/tests/test_dataset.py @@ -4548,18 +4548,18 @@ def test_gradient_datetime(dask): da = da.chunk({'x': 4}) # along x - actual = da.differentiate('x', edge_order=1, time_unit='D') + actual = da.differentiate('x', edge_order=1, datetime_unit='D') expected_x = xr.DataArray( npcompat.gradient( - da, utils.to_numeric(da['x'], time_unit='D'), + da, utils.to_numeric(da['x'], datetime_unit='D'), axis=0, edge_order=1), dims=da.dims, coords=da.coords) assert_equal(expected_x, actual) - actual2 = da.differentiate('x', edge_order=1, time_unit='h') + actual2 = da.differentiate('x', edge_order=1, datetime_unit='h') assert np.allclose(actual, actual2 * 24) # for datetime variable - actual = da['x'].differentiate('x', edge_order=1, time_unit='D') + actual = da['x'].differentiate('x', edge_order=1, datetime_unit='D') assert np.allclose(actual, 1.0) # with different date unit From c31539ec315998d2b699617a893eb10c2e1e1adc Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Thu, 20 Sep 2018 06:53:16 +0900 Subject: [PATCH 23/24] Update via comments --- doc/computation.rst | 4 ++-- doc/whats-new.rst | 4 ++-- xarray/core/dataset.py | 3 +-- 3 files changed, 5 insertions(+), 6 deletions(-) diff --git a/doc/computation.rst b/doc/computation.rst index 125c585edb0..6e484f95801 100644 --- a/doc/computation.rst +++ b/doc/computation.rst @@ -205,7 +205,7 @@ Computation using Coordinates Xarray objects have some handy methods for the computation with their coordinates. :py:meth:`~xarray.DataArray.differentiate` computes derivatives by -finite central differences using their coordinates, +central finite differences using their coordinates, .. ipython:: python a = xr.DataArray([0, 1, 2, 3], dims=['x'], coords=[0.1, 0.11, 0.2, 0.3]) @@ -215,7 +215,7 @@ finite central differences using their coordinates, This method can be used also for multidimensional arrays, .. ipython:: python - a = xr.DataArray(np.arange(8).reshape(4, 2), dims=['x'], + a = xr.DataArray(np.arange(8).reshape(4, 2), dims=['x', 'y'], coords=[0.1, 0.11, 0.2, 0.3]) a.differentiate('x') diff --git a/doc/whats-new.rst b/doc/whats-new.rst index 7365d312353..82d315bdfc8 100644 --- a/doc/whats-new.rst +++ b/doc/whats-new.rst @@ -36,8 +36,8 @@ Documentation Enhancements ~~~~~~~~~~~~ -- :py:func:`~xarray.differentiate`, :py:meth:`~xarray.DataArray.differentiate`, - and :py:meth:`~xarray.Dataset.differentiate` are newly added. +- :py:meth:`~xarray.DataArray.differentiate` and + :py:meth:`~xarray.Dataset.differentiate` are newly added. (:issue:`1332`) By `Keisuke Fujii `_. - Default colormap for sequential and divergent data can now be set via diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py index 5b9044a4a28..d612f96bec9 100644 --- a/xarray/core/dataset.py +++ b/xarray/core/dataset.py @@ -3709,8 +3709,7 @@ def differentiate(self, coord, edge_order=1, datetime_unit=None): for k, v in self.variables.items(): if (k in self.data_vars and dim in v.dims and k not in self.coords): - if v.dtype.kind in 'mM': - v = to_numeric(v, datetime_unit=datetime_unit) + v = to_numeric(v, datetime_unit=datetime_unit) grad = duck_array_ops.gradient( v.data, coord_data, edge_order=edge_order, axis=v.get_axis_num(dim)) From 528bcab00920a40a49643d412ea0d9c8a2d2102c Mon Sep 17 00:00:00 2001 From: fujiisoup Date: Thu, 20 Sep 2018 16:08:30 +0900 Subject: [PATCH 24/24] Update docs. --- doc/computation.rst | 3 +++ xarray/core/dataarray.py | 4 ++++ xarray/core/dataset.py | 4 ++++ 3 files changed, 11 insertions(+) diff --git a/doc/computation.rst b/doc/computation.rst index 6e484f95801..67cda6f2191 100644 --- a/doc/computation.rst +++ b/doc/computation.rst @@ -219,6 +219,9 @@ This method can be used also for multidimensional arrays, coords=[0.1, 0.11, 0.2, 0.3]) a.differentiate('x') +.. note:: + This method is limited to simple cartesian geometry. Differentiation along + multidimensional coordinate is not supported. .. _compute.broadcasting: diff --git a/xarray/core/dataarray.py b/xarray/core/dataarray.py index 8f95912704e..03ea496658f 100644 --- a/xarray/core/dataarray.py +++ b/xarray/core/dataarray.py @@ -2296,6 +2296,10 @@ def differentiate(self, coord, edge_order=1, datetime_unit=None): """ Differentiate the array with the second order accurate central differences. + .. note:: + This feature is limited to simple cartesian geometry, i.e. coord + must be one dimensional. + Parameters ---------- coord: str diff --git a/xarray/core/dataset.py b/xarray/core/dataset.py index d612f96bec9..131345693df 100644 --- a/xarray/core/dataset.py +++ b/xarray/core/dataset.py @@ -3670,6 +3670,10 @@ def differentiate(self, coord, edge_order=1, datetime_unit=None): """ Differentiate with the second order accurate central differences. + .. note:: + This feature is limited to simple cartesian geometry, i.e. coord + must be one dimensional. + Parameters ---------- coord: str