FYI: Xarray's .sum() will convert nan to 0 if all values are nan in a slice, which can affect comparison with CDAT summing results

[tomvothecoder]

Overview

I've ran into this a few times while migrating production code for e3sm_diags and e3sm_to_cmip from CDAT to Xarray/xCDAT

I compare two netCDF, one from Xarray and one from CDAT (post-summing) using np.testing.assert_all_close(). However, I get x and y nan location mismatch because Xarray replaces nan with 0s while CDAT maintains nan (if all values are nan for a slice).

Example:

# Source: https://github.com/E3SM-Project/e3sm_diags/issues/790

Comparing:
/global/cfs/projectdirs/e3sm/e3sm_diags_cdat_test/debug-654-zonal_mean_xy/zonal_mean_xy/Cloud ISCCP/ISCCPCOSP-CLDTOT_TAU1.3_9.4_ISCCP-JJA-global_ref.nc 
 /global/cfs/projectdirs/e3sm/e3sm_diags_cdat_test/main/zonal_mean_xy/Cloud ISCCP/ISCCPCOSP-CLDTOT_TAU1.3_9.4_ISCCP-JJA-global_ref.nc
var_key CLDTOT_TAU1.3_9.4_ISCCP

Not equal to tolerance rtol=1e-05, atol=0

x and y nan location mismatch:
 x: array([ 0.      ,  0.      ,  0.      ,  0.      ,  0.      ,  0.      ,
        0.      ,  0.      ,  0.      , 47.524376, 49.217592, 40.503558,
       40.208695, 42.368686, 41.70946 , 44.171651, 44.821778, 45.508399,...
 y: array([      nan,       nan,       nan,       nan,       nan,       nan,
             nan,       nan,       nan, 48.193733, 49.217592, 40.503558,
       40.208696, 42.368686, 41.70946 , 44.17165 , 44.821777, 45.5084  ,...

Reason why Xarray's .sum() replaces nan with 0

Yes this is by design and matches numpy IIRC. See #2889. I think you have to manually mask the output using .count to get what you want.

Please reopen if you have more questions.
Source: pydata/xarray#5693

Workaround (if you need to compare CDAT vs. Xarray/xCDAT sum results)

1. Use skipna and min_count arg, for example .sum(..., skipna=True, min_count=1).
   • https://docs.xarray.dev/en/stable/generated/xarray.DataArray.sum.html

   • skipna : bool or None, optional
        If True, skip missing values (as marked by NaN). By default, only
        skips missing values for float dtypes; other dtypes either do not
        have a sentinel missing value (int) or ``skipna=True`` has not been
        implemented (object, datetime64 or timedelta64).
     min_count : int or None, optional
        The required number of valid values to perform the operation. If
        fewer than min_count non-NA values are present the result will be
        NA. Only used if skipna is set to True or defaults to True for the
        array's dtype. Changed in version 0.17.0: if specified on an integer
        array and skipna=True, the result will be a float array. 
     

2. In your Xarray/xCDAT code, use np.sum() directly on the .values of the xr.DataArray and rebuild the xr.DataArray object with those results (example).

Related posts:

• Dataset.sum() returns 0 for nan values pydata/xarray#5693
• nansum vs nanmean for all-nan vectors pydata/xarray#2889
• https://numpy.org/doc/stable/reference/generated/numpy.nansum.html
• https://stackoverflow.com/questions/77001292/dealing-with-nan-values-in-xarray-when-using-sum-for-a-basin-netcdf-file

[pochedls]

Just to make sure I understand: If an array has NaN values, xarray will return the sum of all non-NaN values and CDMS/numpy will return NaN. Right?

[lee1043]

@pochedls I think your understanding is correct, according to the xarray's sum.

[tomvothecoder]

Xarray's .sum() matches numpy's behavior when using np.nansum(), which treats nan as zero.

Thanks @mathause I also think the current behavior is not perfect but the best.

I would expect both to return np.nan

I expect that np.nansum(ds) is equivalent to np.sum(not nan values) and thus should be 0, while np.mean should be NaN as @dcherian pointed out.

To me, the future average function would also return np.nan for all nan slices.

Source - pydata/xarray#2889

Return the sum of array elements over a given axis treating Not a Numbers (NaNs) as zero.

In NumPy versions <= 1.9.0 Nan is returned for slices that are all-NaN or empty. In later versions zero is returned.

Source - https://numpy.org/doc/stable/reference/generated/numpy.nansum.html

[tomvothecoder]

If all values for a subset are np.nan, then the sum will be 0. This causes issues when comparing against CDAT, which uses mask and maintains np.nan.

[jypeter]

Yet another example of what happens when you rely on NaNs, rather than on masks ?

Though I confess that the documentation of ma.MaskedArray.sum is not very detailed. You need to know that the operations will only take the non-masked elements into account, but this is also what you expect

[tomvothecoder]

Refer to my comment above for more info on why this behavior happens (it's not wrong).