code: def interpolate_xarray_linear(xpoints, ypoints, values, shape, chunks=CHUNK_SIZE):
      
      from scipy.interpolate.interpnd import (LinearNDInterpolator,
                                              _ndim_coords_from_arrays)
  
      if isinstance(chunks, (list, tuple)):
          vchunks, hchunks = chunks
      else:
          vchunks, hchunks = chunks, chunks
  
      points = _ndim_coords_from_arrays(np.vstack((np.asarray(ypoints),
                                                   np.asarray(xpoints))).T)
  
      interpolator = LinearNDInterpolator(points, values)
  
      grid_x, grid_y = da.meshgrid(da.arange(shape[1], chunks=hchunks),
                                   da.arange(shape[0], chunks=vchunks))
  
      # workaround for non-thread-safe first call of the interpolator:
      interpolator((0, 0))
      res = da.map_blocks(intp, grid_x, grid_y, interpolator=interpolator)
  
      return DataArray(res, dims=('y', 'x'))
docs: Interpolate linearly, generating a dask array.