# Possible bug(?) in xarray.groupby('time.month') operations **URL:** \ **Category:** Technical\ **Tags:** help, xarray, jupyterlab\ **Created:** [28 January 2025 04:50 UTC](https://forum.access-hive.org.au/t/possible-bug-in-xarray-groupby-time-month-operations/4101 "2025-01-28T04:50:43Z")\ **Posts on this page:** 6\ **Page:** 2
**Author:** ![Thomas-Moore](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/thomas-moore/32/329_2.png) [@Thomas-Moore](https://forum.access-hive.org.au/u/Thomas-Moore)\ **Post date:** [31 January 2025 02:26 UTC](https://forum.access-hive.org.au/t/possible-bug-in-xarray-groupby-time-month-operations/4101/21 "2025-01-31T02:26:53Z")
Thanks for these efforts @jemmajeffree et al. FYI for all - if you want a good, basic overview of Xarray’s “groupby” from Deepak spend 37 minutes on this \> [https://youtu.be/92-QU37W9WI?si=a1XgLxqsjygjPyfu](https://youtu.be/92-QU37W9WI?si=a1XgLxqsjygjPyfu) ---
**Author:** ![navidcy](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/navidcy/32/101_2.png) [@navidcy](https://forum.access-hive.org.au/u/navidcy)\ **Post date:** [4 February 2025 01:44 UTC](https://forum.access-hive.org.au/t/possible-bug-in-xarray-groupby-time-month-operations/4101/22 "2025-02-04T01:44:15Z")
I think this needs to be posted at xarray repo? I admit I don’t know what flox is… But the example above with groupby giving NaNs seems very counter-intuitive to me (as did the initial post by @hrsdawson)!! edit: I actually missed reading the “I’m planning on raising the bug with the xarray…” bit by @jemmajeffree; nice! ---
**Author:** ![navidcy](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/navidcy/32/101_2.png) [@navidcy](https://forum.access-hive.org.au/u/navidcy)\ **Post date:** [5 February 2025 16:35 UTC](https://forum.access-hive.org.au/t/possible-bug-in-xarray-groupby-time-month-operations/4101/23 "2025-02-05T16:35:53Z")
> [@jemmajeffree](#): > > _My best guess as to what’s happening is that calling std after groupby on a dask array is lower precision than a numpy array, which is producing noisy variances, and occasionally negative variances that produce NaNs when square rooted_ Well **if indeed** somehow `groupby` results in the computations following it being done at float32 instead of float64 that might explain things! For datasets where the variation is close to the float32 presicion then adding/subtracting/etc gives nonsense (including variances that are negative) and thus NaN when `sqrt` is taken. ---
**Author:** ![Benoit](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/benoit/32/3115_2.png) [@Benoit](https://forum.access-hive.org.au/u/Benoit)\ **Post date:** [13 February 2025 23:40 UTC](https://forum.access-hive.org.au/t/possible-bug-in-xarray-groupby-time-month-operations/4101/24 "2025-02-13T23:40:11Z")
I’m just lurking but I am curious (and a bit anxious). Was an issue posted on GitHub about this in the end? Could it be related to issues such as [bottleneck : Wrong mean for float32 array · Issue #1346 · pydata/xarray · GitHub](https://github.com/pydata/xarray/issues/1346)? If using Dask toggles a precision switch that silently returns wrong outputs this is a pretty big deal, right? ---
**Author:** ![jemmajeffree](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/jemmajeffree/32/228_2.png) [@jemmajeffree](https://forum.access-hive.org.au/u/jemmajeffree)\ **Post date:** [14 February 2025 06:20 UTC](https://forum.access-hive.org.au/t/possible-bug-in-xarray-groupby-time-month-operations/4101/25 "2025-02-14T06:20:49Z")
Hi Benoit and others. I haven’t raised the issue yet, but I’ve identified the lines causing the problem and am partway through implementing a fix. It’s not toggling a precision switch, it’s just a different numerical implementation of standard deviation that gets noisy when the standard deviation is tiny compared to the mean (i.e deep ocean salinity). Because I had to finish something for a deadline yesterday, I only got a chance to dive into the code this morning. The discrepancy comes down to flox using [these lines](https://github.com/xarray-contrib/flox/blob/ca576812e78b3978421eace6e9dde5a76729ebcc/flox/aggregate_npg.py#L112) with numpy/loaded data and [these lines](https://github.com/xarray-contrib/flox/blob/ca576812e78b3978421eace6e9dde5a76729ebcc/flox/aggregations.py#L379) when using dask arrays and lazy data; the numpy version starts by subtracting an offset from the array so there isn’t such a huge difference between the magnitude of mean and standard deviation. I’m in the process of implementing the same step into the dask implementation of flox. It’s working for the single-dimension case but I still need to generalise to work with any number of dimensions. Once I’ve got this done (hopefully Monday) I’ll post the issue and fix on github ---
**Author:** ![jemmajeffree](https://sea2.discourse-cdn.com/flex020/user_avatar/forum.access-hive.org.au/jemmajeffree/32/228_2.png) [@jemmajeffree](https://forum.access-hive.org.au/u/jemmajeffree)\ **Post date:** [17 February 2025 01:29 UTC](https://forum.access-hive.org.au/t/possible-bug-in-xarray-groupby-time-month-operations/4101/26 "2025-02-17T01:29:04Z")
I’ve now raised the issue on the flox repo: > > > Hi, > I've noticed that in a few rare situations, groupby and flox can return quit…e noisy standard deviations. In situations where the mean of an array is much larger than the standard deviation (such as deep ocean salinity, raised \[here\]()), flox returns noisier values on dask arrays than on loaded numpy arrays. In extreme situations, the standard deviation of a dask array can contain NaNs from square-rooting negative variances. > > I'm guessing it's the same idea as #386, in which case @dcherian has thought about this for much longer than I have. I've done a little bit of looking through the code, and could easily have missed something about how this works with neighbouring functions, but my thoughts on the potential problem and how it might be addressed are below. > > Minimal complete verifiable example: > \`\`\`python > import numpy as np > import xarray as xr > > l =12000 > np.random.seed(1) > test\_data = xr.DataArray(np.random.uniform(0,1,l)/100+1000000,dims=('time',) # huge mean with relatively small variability > ).assign\_coords({'month':xr.DataArray(np.arange(l)%12,dims=('time',))}) > > \# with numpy arrays returns reasonable and consistent values > test\_data.groupby('month').std('time') > \# array(\[0.00283648, 0.00281895, 0.00287791, 0.00287652, 0.00287337, > \# 0.00287037, 0.00289802, 0.00289441, 0.00285839, 0.00296478, > \# 0.00284787, 0.00292089\]) > > \# using lazy computation/dask > dask\_test\_data = test\_data.chunk({'time':100}) > dask\_test\_data.groupby('month').std('time').load() > \# array(\[0.01118034, 0.01118034, 0.01118034, 0.01581139, 0. , > \# 0.01581139, 0.01118034, 0.01118034, 0.01118034, nan, > \# nan, 0. \]) > \`\`\` > > A functional workaround is to subtract the mean before calculating standard deviation: > \`\`\`python > (dask\_test\_data.groupby('month')-dask\_test\_data.groupby('month').mean('time')).groupby('month').std('time').load() > \`\`\` > > My understanding is that the distinction comes from aggregate\_npg.py improving precision by \[subtracting the first non-nan element of the array\](https://github.com/xarray-contrib/flox/blob/ca576812e78b3978421eace6e9dde5a76729ebcc/flox/aggregate\_npg.py#L112), a preprocessing step \[skipped by aggregations.py\](https://github.com/xarray-contrib/flox/blob/ca576812e78b3978421eace6e9dde5a76729ebcc/flox/aggregations.py#L379). This solution is probably not quite as stable as subtracting the mean, but the first element should be really close to the mean if the standard deviation is small, and it might be faster. > > I’d suggest that to improve precision and match the numpy engine behaviour in aggregations\_npg.py, the flox engine implementation for dask arrays of nanstd,nanvar,std,var could have a preprocessor that looks something like this: > > \`\`\`python > def var\_std\_preprocess(array, axis): # Not sure of naming conventions, sorry > """Subtracts first value of array from whole array, > to improve numerical precision of nanstd, nanvar, std, var > > Adapted from from argreduce\_preprocess and \_var\_std\_wrapper in aggregate\_npg.py > """ > import dask.array # Copied from argreduce\_preprocess, but maybe these shouldn’t be within the function? > import numpy as np # For either this function or argreduce\_preprocess? > > # NEXT LINE IS PSEUDOCODE; I’m not entirely sure how to apply it lazily > # If it doesn’t cost anything speed wise, then probably better to use mean. Happy to run some time tests on either > first\_elements = nanfirst(array,axis) > > def subtract\_first(array\_, first\_elements\_): > return array\_-first\_elements\_ > > return dask.array.map\_blocks( > subtract\_first, > array, > first, > dtype=array.dtype, > meta=array.\_meta, > name="groupby-var\_std-preprocess", > ) > \`\`\` > > and is included in the Aggregations definition like so: > \`\`\`python > nanstd = Aggregation( > "nanstd", > preprocess=var\_std\_preprocess, #UPDATED LINE > chunk=("nansum\_of\_squares", "nansum", "nanlen"), > combine=("sum", "sum", "sum"), > finalize=\_std\_finalize, > fill\_value=0, > final\_fill\_value=np.nan, > dtypes=(None, None, np.intp), > final\_dtype=np.floating, > ) > \`\`\` > > It seems to work if \`first\_elements\` is naively \`array\[0\]\` in the one-dimensional, no-nans case, but I’m not sure how to generalise it and apply nanfirst without the usual layers/wrappers. (aggregate\_npg.py uses \`first = \_get\_aggregate(engine).aggregate(group\_idx, array, func="nanfirst", axis=axis)\`, but I don't think this syntax translates to the flox/dask implementation) . If you can give me a few tips or examples to work from, then I’m happy to try implement this behaviour. > > Happy also to discuss alternatives, or to provide a pull request if that's easier to work with. > > This is also my first time reading through flox code in detail (it's really nicely written and documented, by the way, was lovely to read), and one of the first times I’ve interacted with public github repos, so I’d appreciate any feedback or corrections on what's useful to provide when describing issues. I didn’t quite manage to generalise to multiple dimensions, hoping someone will give me a hand with that. For the moment, subtract the mean before taking standard deviations with groupby, if the mean is big compared to the expected standard deviation [Previous page](https://forum.access-hive.org.au/t/possible-bug-in-xarray-groupby-time-month-operations/4101.md?page=1)