[Python] precision in numpy versus xarray

In a recent analysis, I stumbled across some surprising behavior of my data. I (think I've) resolved it by converting my data to doubles but I wanted to bring this up here to find out if this is behaving as it should be and/or to satisfy my curiosity to why this is happening.

I've simplified things to the example below. Basically there is quite a difference in the result when using float32 between numpy and xarray, where I thought that under the hood the operation that is happening should be identical (for my actual case, this becomes way more serious):

import numpy as np
import xarray as xr

samples = 10_000_000
data = np.random.RandomState(123).chisquare(2, size=samples)
data = data.astype(np.float32)
da = xr.DataArray(data)

print('numpy: {:e}'.format((data - data.mean()).mean()))
print('xarray: {:e}'.format((da - da.mean()).mean().item()))


numpy: -8.652336e-07
xarray: 2.263663e-02


This is how I solved my issues:

data = data.astype(np.double)
da = xr.DataArray(data)

print('numpy: {:e}'.format((data - data.mean()).mean()))
print('xarray: {:e}'.format((da - da.mean()).mean().item()))


numpy: -1.637659e-15
xarray: 6.976170e-15

Continue reading...