Problem with boundary values in .core.grids._horizontal_gradient

[lesommer]

numpy.gradient compute the spatial derivatives with a second order scheme except near boundaries where the derivatives are evaluated with a first order scheme. oocgcm.core.grids._horizontal_gradient wraps numpy.gradient for xarray.DataArray with dask.

in order for the result not to depend on the chunk size, we need to set depth and use dask.array.map_overlap.

but it is not clear how to specify that the calculation should be first order near the boundary when the chunk is close to the boundary. numpy.gradient sees the ghost cell and therefore applies the second order derivative...

the current implementation explicitely sets a np.nan boundary condition to avoid errors.

this could be a question for dask mailing list.

[kuchaale]

I think your current version would not work anyway if you don't specify a dtype explicitly, see my issue here.

[lesommer]

@kuchaale : thanks for your comment and your suggestion. But I am not sure to get what you mean by "would not work" here. The present version of the code actually produce a result, but does not compute the gradient close to the boundary (thanks to the specification of np.nan).

[kuchaale]

@lesommer : I understand your issue with the specification of np.nan. However, when I run the code below (without dtype specification):

`import numpy as np
import xarray as xr

def is_numpy(array):
test = bool( isinstance(array,np.ndarray)
+ isinstance(array,np.ma.masked_array) )
return test

def _horizontal_gradient(scalararray):
data = scalararray.data
coords = scalararray.coords
dims = scalararray.dims
chunks = scalararray.chunks
if is_numpy(data):
da_dj,da_di = np.gradient(data)
else:
x_derivative = lambda arr:np.gradient(arr,axis=-1) # req. numpy > 1.11
y_derivative = lambda arr:np.gradient(arr,axis=-2) # req. numpy > 1.11
gx = data.map_overlap(x_derivative,depth=(0,1),boundary={1: np.nan})#, dtype=data.dtype)
gy = data.map_overlap(y_derivative,depth=(1,0),boundary={0: np.nan})#, dtype=data.dtype)
da_di = xr.DataArray(gx,coords,dims)
da_dj = xr.DataArray(gy,coords,dims)
return da_dj,da_di

ds = xr.open_dataset('../data/heigh_2017011715Z_1000hPa.nc')
ds = ds.chunk(chunks = 40)
ds`

I get this error:
ValueError: dtypeinference failed inmap_blocks`.

Original error is below:

ValueError('Shape of array too small to calculate a numerical gradient, at least two elements are required.',)

Traceback:

File "c:\python27\lib\site-packages\dask\array\core.py", line 457, in apply_infer_dtype
o = func(*args, **kwargs)
File "", line 26, in
x_derivative = lambda arr:np.gradient(arr,axis=-1) # req. numpy > 1.11
File "c:\python27\lib\site-packages\numpy\lib\function_base.py", line 1642, in gradient
"Shape of array too small to calculate a numerical gradient, "
`
This issue is induced by the non-specification of dtype.

[redouanelg]

Hi @kuchaale, your solution worked for me, please open a PR or allow me to do it