Some tips if you're comparing Xarray/xCDAT objects against cdms2 Masked Arrays

[tomvothecoder]

1. Make sure representation of missing values are the same between cdms2 and Xarray

• cdms2 variables use numpy masked arrays. Missing values are represented by a fill_value set in the NetCDF file (e.g., fill_value=1e-20) or manually set by the user.
• Xarray uses numpy arrays and represents missing values with np.nan (docs).

The problem: If missing values are represented differently, then there can be significant differences in floating point
comparisons.

Solutions:

1. Replace missing values in cdms2 variable objects with np.nan to align with Xarray's behavior with missing values.

   # t_var is a cdms2 TransientVariable and `.filled()` is a `numpy.ma` method.
   # https://numpy.org/doc/stable/reference/generated/numpy.ma.filled.html
   t_var = t_var.filled(np.nan)

2. Replace Xarray and cdms2 variable object missing values with the same fill value (e.g., 1e-20).

   # `.fillna()` is an xarray DataArray method
   # https://docs.xarray.dev/en/stable/generated/xarray.DataArray.fillna.html
   dataarray = dataarray.fillna(1e-20)
   
   # `.filled()` is a numpy MaskedArray method
   # https://numpy.org/doc/stable/reference/generated/numpy.ma.filled.html
   t_var = t_var.filled(1e-20)

2. Compare the underlying numpy arrays found in cdms2 and Xarray data structures

cdms2 and Xarray objects are both interoperable with np.testing.assert_allclose(). However, I'm not 100% sure how cdms2 masked arrays behave when passed to np.testing.assert_allclose().

Extracting the underlying numpy arrays ensures that cdms2 and/or Xarray does not introduce potential unwanted side-effects when comparing values between packages.

Example:

tvar = cdms2.TransientVariable(...)
da = xr.DataArray(...)

# .data is the underlying numpy array found in both data structures
np.testing.assert_allclose(tvar.data, da.data)

3. Use np.testing.assert_allclose() to compare max absolute difference (floating point) and max relative difference (%).

• Comparing absolute differences makes sense if the numbers are not small floating points and have consistent digit places (e.g., [100.01, 95.02, ...])
• Comparing relative difference makes sense if the numbers in an array are small floating points or vary significant in digit places (e.g., [0.000123, 0.000959, ...]). This can apply in scenarios such as comparing climatological anomalies.

Related resources:

• https://floating-point-gui.de/errors/comparison/
• https://bitbashing.io/comparing-floats.html
• https://stackoverflow.com/questions/65909842/what-is-rtol-for-in-numpys-allclose-function

[tomvothecoder]

Hey @gleckler1, if I recall correctly, you mentioned running into issues with large differences when comparing CDAT vs. Xarray/xCDAT results. These tips might be helpful.

[jypeter]

Thanks @tomvothecoder !

For people like me who don't like np.nan, we can probably do the same stuff by masking the nan values with np.ma.masked_invalid. Fortunately, testing.assert_allclose seems to work on both numpy array with nans or masks

And, if A is a cdms2 masked variable, you can get the pure np.ma array part (that is, no axes, metadata, etc...) with A.asma()

We have to b careful not to mix masked arrays with np arrays with nans, or it will silently get ugly :-(

>>> a_np = np.identity(3)
>>> a_np[:, 0] = np.nan
>>> a_np
array([[nan,  0.,  0.],
       [nan,  1.,  0.],
       [nan,  0.,  1.]])
>>> a_ma = np.ma.masked_invalid(a_np, copy=True)
>>> a_ma
masked_array(
  data=[[--, 0.0, 0.0],
        [--, 1.0, 0.0],
        [--, 0.0, 1.0]],
  mask=[[ True, False, False],
        [ True, False, False],
        [ True, False, False]],
  fill_value=1e+20)
>>> np.testing.assert_allclose(a_ma, a_np)
>>> 
>>> b_np = np.identity(3)
>>> b_np[0, :] = np.nan
>>> b_np
array([[nan, nan, nan],
       [ 0.,  1.,  0.],
       [ 0.,  0.,  1.]])
>>> a_ma + b_np
masked_array(
  data=[[--, nan, nan],
        [--, 2.0, 0.0],
        [--, 0.0, 2.0]],
  mask=[[ True, False, False],
        [ True, False, False],
        [ True, False, False]],
  fill_value=1e+20)

>>> a_ma.count()
6
>>> (a_ma + b_np).count()
6

[tomvothecoder]

Thank you @jypeter, I haven't tried np.ma.masked_invalid but that looks like a useful alternative if you don't want to work with np.nan.

Good point about not mixing masked_invalid elements with np.nan elements (np.nan won't be considered missing).