Memory issue merging NetCDF files using xarray.open_mfdataset and to_netcdf

What happened?

I have 5 NetCDF files (1 GiB each). They have 4 dimensions: time, depth, lat, lon. All the files have exactly the same depth, lat, lon. The time axis have the same interval and there are no gaps on this axis for all the 5 files (and there is continuity in the axis between files).

All I am doing is merging the files along the time-axis and saving it to a new NetCDF file.

Running the script, I allocated 185 GiB of memory (the maximum in my cluster).

The program runs until the to_netcdf() function is called. I get an error stating there is not enough memory.

What did you expect to happen?

As the 5 files are 1 GiB each, and I allocated 185 GiB (far more than 5² GiB), I expected the program to run and not require more than the allocated memory (after all, I gave 37 times the combined size of the files).

Minimal Complete Verifiable Example

path = './data/data_*.nc' # files are: data_1.nc data_2.nc data_3.nc data_4.nc data_5.nc
data = xr.open_mfdataset(path)

data = data.load() # uses 5 GiB - tested with a memory profiler

data.to_netcdf('./output/combined.nc') #

MVCE confirmation

• Minimal example — the example is as focused as reasonably possible to demonstrate the underlying issue in xarray.
• Complete example — the example is self-contained, including all data and the text of any traceback.
• Verifiable example — the example copy & pastes into an IPython prompt or Binder notebook, returning the result.
• New issue — a search of GitHub Issues suggests this is not a duplicate.

Relevant log output

Traceback (most recent call last):
  File "/users/me/code/par2.py", line 78, in <module>
    preprocess_data(year, month)
  File "/users/me/code/par2.py", line 69, in preprocess_data
    data.to_netcdf(path=outpath)
  File "/CSC_CONTAINER/miniconda/envs/env1/lib/python3.10/site-packages/xarray/core/dataset.py", line 1882, in to_netcdf
    return to_netcdf(  # type: ignore  # mypy cannot resolve the overloads:(
  File "/CSC_CONTAINER/miniconda/envs/env1/lib/python3.10/site-packages/xarray/backends/api.py", line 1210, in to_netcdf
    dump_to_store(
  File "/CSC_CONTAINER/miniconda/envs/env1/lib/python3.10/site-packages/xarray/backends/api.py", line 1257, in dump_to_store
    store.store(variables, attrs, check_encoding, writer, unlimited_dims=unlimited_dims)
  File "/CSC_CONTAINER/miniconda/envs/env1/lib/python3.10/site-packages/xarray/backends/common.py", line 263, in store
    variables, attributes = self.encode(variables, attributes)
  File "/CSC_CONTAINER/miniconda/envs/env1/lib/python3.10/site-packages/xarray/backends/common.py", line 352, in encode
    variables, attributes = cf_encoder(variables, attributes)
  File "/CSC_CONTAINER/miniconda/envs/env1/lib/python3.10/site-packages/xarray/conventions.py", line 864, in cf_encoder
    new_vars = {k: encode_cf_variable(v, name=k) for k, v in variables.items()}
  File "/CSC_CONTAINER/miniconda/envs/env1/lib/python3.10/site-packages/xarray/conventions.py", line 864, in <dictcomp>
    new_vars = {k: encode_cf_variable(v, name=k) for k, v in variables.items()}
  File "/CSC_CONTAINER/miniconda/envs/env1/lib/python3.10/site-packages/xarray/conventions.py", line 273, in encode_cf_variable
    var = coder.encode(var, name=name)
  File "/CSC_CONTAINER/miniconda/envs/env1/lib/python3.10/site-packages/xarray/coding/variables.py", line 170, in encode
    data = duck_array_ops.fillna(data, fill_value)
  File "/CSC_CONTAINER/miniconda/envs/env1/lib/python3.10/site-packages/xarray/core/duck_array_ops.py", line 283, in fillna
    return where(notnull(data), data, other)
  File "/CSC_CONTAINER/miniconda/envs/env1/lib/python3.10/site-packages/xarray/core/duck_array_ops.py", line 270, in where
    return _where(condition, *as_shared_dtype([x, y]))
  File "<__array_function__ internals>", line 180, in where
numpy.core._exceptions._ArrayMemoryError: Unable to allocate 43.6 GiB for an array with shape (280, 200, 277, 754) and data type float32

Anything else we need to know?

I allocated 185 GiB for this job, from my understanding, this means that merging 5 datasets with 1 GiB each requires more than 185 GiB memory. It sounds like a memory leak to me.

I am not the only one with this issue, cf: https://github.com/pydata/xarray/discussions/4890

Environment

Start with the MVCE and traceback through xarray.Dataset.to_netcdf in xarray/core/dataset.py, then follow the calls in backends/api.py, conventions.py, and coding/variables.py. Reproduce the allocation failure and determine the expected memory behavior; done should include a documented resolution or a focused regression test for the observed path.