diff() failing on lazy-loaded data only

[lluritu]

What happened?

Hello,
I am reading some GRIB2 data from ECMWF MARS archive and cannot diff("step") unless I load the data previously. I am seeing the problem only after upgrading from xarray v2025.3.1 to v2025.10.1

xr.open_dataset("/path/to/tp_2025052500.grb").load().diff("step") works well.
xr.open_dataset("/path/to/tp_2025052500.grb").diff("step") fails with alignment problems.

Anyone any clues what's going on?
Thanks
Llorenç

[lluritu]

An addition:
xr.open_dataset("/path/to/tp_2025052500.grb", chunks='auto').diff("step") also works well.
So the problem is narrowed down to xarray Lazy Indexing classes (used when chunks=None).

Llorenç

[keewis]

can you share a link to the file (if it is public), so we can try to reproduce this? Also, could you please share the traceback?

[lluritu]

You can get some sample data from ECMWF open data archive for two different steps using wget and concat them together:

wget https://data.ecmwf.int/forecasts/20251121/00z/ifs/0p25/oper/20251121000000-0h-oper-fc.grib2
wget https://data.ecmwf.int/forecasts/20251121/00z/ifs/0p25/oper/20251121000000-12h-oper-fc.grib2
cat 20251121000000-0h-oper-fc.grib2 20251121000000-12h-oper-fc.grib2 > data.grib2

The data.grib2 will contain two steps.
Note: this particular files will only be available for a few days, but some recent ones are always there.
Those files contain many parameters, might be useful to restrict to a single parameter to simplify. In my case I was testing total precipitation, but I don't think the particular parameter matters.

[lluritu]

Traceback:

---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
Cell In[200], line 1
----> 1 xr.open_dataset("[/path/to/cat.grib2]").diff("step")

File [/usr/local/apps/python3/3.12.11-01/lib/python3.12/site-packages/xarray/core/dataset.py:7835](https://jupyterhub.ecmwf.int/usr/local/apps/python3/3.12.11-01/lib/python3.12/site-packages/xarray/core/dataset.py#line=7834), in Dataset.diff(self, dim, n, label)
   7833 elif dim in var.dims:
   7834     if name in self.data_vars:
-> 7835         variables[name] = var.isel(slice_end) - var.isel(slice_start)
   7836     else:
   7837         variables[name] = var.isel(slice_new)

File [/usr/local/apps/python3/3.12.11-01/lib/python3.12/site-packages/xarray/core/_typed_ops.py:920](https://jupyterhub.ecmwf.int/usr/local/apps/python3/3.12.11-01/lib/python3.12/site-packages/xarray/core/_typed_ops.py#line=919), in VariableOpsMixin.__sub__(self, other)
    919 def __sub__(self, other: VarCompatible) -> Self | T_DA | Dataset | DataTree:
--> 920     return self._binary_op(other, operator.sub)

File [/usr/local/apps/python3/3.12.11-01/lib/python3.12/site-packages/xarray/core/variable.py:2437](https://jupyterhub.ecmwf.int/usr/local/apps/python3/3.12.11-01/lib/python3.12/site-packages/xarray/core/variable.py#line=2436), in Variable._binary_op(self, other, f, reflexive)
   2435     other_data, self_data, dims = _broadcast_compat_data(other, self)
   2436 else:
-> 2437     self_data, other_data, dims = _broadcast_compat_data(self, other)
   2438 keep_attrs = _get_keep_attrs(default=False)
   2439 attrs = self._attrs if keep_attrs else None

File [/usr/local/apps/python3/3.12.11-01/lib/python3.12/site-packages/xarray/core/variable.py:3046](https://jupyterhub.ecmwf.int/usr/local/apps/python3/3.12.11-01/lib/python3.12/site-packages/xarray/core/variable.py#line=3045), in _broadcast_compat_data(self, other)
   3038     raise ValueError(
   3039         "Broadcasting is necessary but automatic broadcasting is disabled via "
   3040         "global option `'arithmetic_broadcast'`. "
   3041         "Use `xr.set_options(arithmetic_broadcast=True)` to enable automatic broadcasting."
   3042     )
   3044 if all(hasattr(other, attr) for attr in ["dims", "data", "shape", "encoding"]):
   3045     # `other` satisfies the necessary Variable API for broadcast_variables
-> 3046     new_self, new_other = _broadcast_compat_variables(self, other)
   3047     self_data = new_self.data
   3048     other_data = new_other.data

File [/usr/local/apps/python3/3.12.11-01/lib/python3.12/site-packages/xarray/core/variable.py:3012](https://jupyterhub.ecmwf.int/usr/local/apps/python3/3.12.11-01/lib/python3.12/site-packages/xarray/core/variable.py#line=3011), in _broadcast_compat_variables(*variables)
   3006 def _broadcast_compat_variables(*variables):
   3007     """Create broadcast compatible variables, with the same dimensions.
   3008 
   3009     Unlike the result of broadcast_variables(), some variables may have
   3010     dimensions of size 1 instead of the size of the broadcast dimension.
   3011     """
-> 3012     dims = tuple(_unified_dims(variables))
   3013     return tuple(var.set_dims(dims) if var.dims != dims else var for var in variables)

File [/usr/local/apps/python3/3.12.11-01/lib/python3.12/site-packages/xarray/core/variable.py:2999](https://jupyterhub.ecmwf.int/usr/local/apps/python3/3.12.11-01/lib/python3.12/site-packages/xarray/core/variable.py#line=2998), in _unified_dims(variables)
   2997             all_dims[d] = s
   2998         elif all_dims[d] != s:
-> 2999             raise ValueError(
   3000                 "operands cannot be broadcast together "
   3001                 f"with mismatched lengths for dimension {d!r}: {(all_dims[d], s)}"
   3002             )
   3003 return all_dims

ValueError: operands cannot be broadcast together with mismatched lengths for dimension 'step': (1, 0)

[dcherian]

Can you show us the repr for xr.open_dataset("[/path/to/cat.grib2]") please?

[lluritu]

This is the repr:

<xarray.Dataset> Size: 1GB
Dimensions:             (step: 2, soilLayer: 4, latitude: 721, longitude: 1440,
                         isobaricInhPa: 13)
Coordinates: (12/13)
  * step                (step) timedelta64[ns] 16B 00:00:00 12:00:00
  * soilLayer           (soilLayer) float64 32B 1.0 2.0 3.0 4.0
  * latitude            (latitude) float64 6kB 90.0 89.75 89.5 ... -89.75 -90.0
  * longitude           (longitude) float64 12kB -180.0 -179.8 ... 179.5 179.8
  * isobaricInhPa       (isobaricInhPa) float64 104B 1e+03 925.0 ... 100.0 50.0
    time                datetime64[ns] 8B ...
    ...                  ...
    entireAtmosphere    float64 8B ...
    heightAboveGround   float64 8B ...
    surface             float64 8B ...
    nominalTop          float64 8B ...
    meanSea             float64 8B ...
    mostUnstableParcel  float64 8B ...
Data variables: (12/39)
    sot                 (step, soilLayer, latitude, longitude) float32 33MB ...
    q                   (step, isobaricInhPa, latitude, longitude) float32 108MB ...
    tcw                 (step, latitude, longitude) float32 8MB ...
    d2m                 (step, latitude, longitude) float32 8MB ...
    ssrd                (step, latitude, longitude) float32 8MB ...
    tcwv                (step, latitude, longitude) float32 8MB ...
    ...                  ...
    svn                 (step, latitude, longitude) float32 8MB ...
    zos                 (step, latitude, longitude) float32 8MB ...
    skt                 (step, latitude, longitude) float32 8MB ...
    tcc                 (step, latitude, longitude) float32 8MB ...
    t2m                 (step, latitude, longitude) float32 8MB ...
    sp                  (step, latitude, longitude) float32 8MB ...
Attributes:
    GRIB_edition:            2
    GRIB_centre:             ecmf
    GRIB_centreDescription:  European Centre for Medium-Range Weather Forecasts
    GRIB_subCentre:          0
    Conventions:             CF-1.7
    institution:             European Centre for Medium-Range Weather Forecasts
    history:                 2025-11-24T08:50 GRIB to CDM+CF via cfgrib-0.9.1...

[keewis]

the problem is this:

xarray/xarray/core/indexing.py

Lines 366 to 368 in 42669f9

	if is_full_slice(old_indexer):
	# shortcut for full slices
	return applied_indexer

(introduced in #10580)

This means that the slices are not normalized anymore if we're indexing the entire dimension, and apparently cfgrib's backend array does not support unnormalized slices? The fix is change the second line to

return normalize_slice(applied_indexer, size)