Finite differencing with xarray and dask.ipynb

Hi @spencerkclark,
I like your approach very much. However, when I tried to produce similar results to ncl example 1 using center_finite_diff_n (see below), I receive different results especially on the boundary points. See below:

t = np.arange(6); x = [30.,33.,39.,36.,41.,37.]; test = xr.DataArray(x, coords=[t], dims=['t'])

Using xdiff

xdiff(test, 't', accuracy=2, method='centered', spacing=2.)
produced:
`<xarray.DataArray (t: 6)>
array([-1. , 2.25, 0.75, 0.5 , 0.25, -2.75])
Coordinates:

• t (t) int64 0 1 2 3 4 5`

Using xgradient

xgradient(test, 't', accuracy=2, spacing=2.).compute()
produced:
`<xarray.DataArray (t: 6)>
array([ 0.75, 5.25, 0.75, 0.5 , 4.5 , -4.25])
Coordinates:

• t (t) int64 0 1 2 3 4 5`

Am I doing something wrong? Thx for your answer,
Ales

Hi @kuchaale,

Sorry I just saw this now (for some reason I wasn't notified by your tagging of me in your post). Many thanks for trying things out!

For right now I would say the boundary behavior of xdiff you experienced is expected, because I've only allowed for the mode='wrap' option, which is enabled by default, and the interior values look good.

But it's clear you've uncovered a bug in xgradient. I think I've fixed things in the updated gist:

In [1]: test = xr.DataArray([30., 33., 39., 36., 41., 37.], dims=['x'])
Out [1]: <xarray.DataArray (x: 6)>
array([-1.  ,  2.25,  0.75,  0.5 ,  0.25, -2.75])
Dimensions without coordinates: x

In [2]: xgradient(test, 'x', spacing=2., accuracy=2)
Out [2]: <xarray.DataArray (x: 6)>
array([ 0.75,  2.25,  0.75,  0.5 ,  0.25, -4.25])
Dimensions without coordinates: x

Note that my result using xgradient is different than NCL's on the boundaries, because I'm doing second-order accurate differencing everywhere. NCL uses second-order centered differencing on the interior and first-order forward or backward differencing on the boundaries.

In the updated gist I also added some logic in xgradient that checks to see if the underlying data is stored in a dask array; only then will it chunk the left and right components. (This way you don't need to call .compute() at the end anymore).