xarray: slow performance with open_mfdataset

In your twitter thread you said

Do any of my xarray/dask folks know why open_mfdataset takes such a significant amount of time compared to looping over a list of files? Each file corresponds to a new time, just wanting to open multiple times at once…

The general reason for this is usually that open_mfdataset performs coordinate compatibility checks when it concatenates the files. It’s useful to actually read the code of open_mfdataset to see how it works.

First, all the files are opened individually https://github.com/pydata/xarray/blob/577d3a75ea8bb25b99f9d31af8da14210cddff78/xarray/backends/api.py#L900-L903

You can recreate this step outside of xarray yourself by doing something like

from glob import glob
datasets = [xr.open_dataset(fname, chunks={}) for fname in glob('*.nc')]

Once each dataset is open, xarray calls out to one of its combine functions. This logic has gotten more complex over the years as different options have been introduced, but the gist is this: https://github.com/pydata/xarray/blob/577d3a75ea8bb25b99f9d31af8da14210cddff78/xarray/backends/api.py#L947-L952

You can reproduce this step outside of xarray, e.g.

ds = xr.concat(datasets, dim='time')

At that point, various checks will kick in to be sure that the coordinates in the different datasets are compatible. Performing these checks requires the data to be read eagerly, which can be a source of slow performance.

Without seeing more details about your files, it’s hard to know exactly where the issue lies. A good place to start is to simply drop all coordinates from your data as a preprocessing step.

def drop_all_coords(ds):
    return ds.reset_coords(drop=True)

xr.open_mfdataset('*.nc', combine='by_coords', preprocess=drop_all_coords)

If you observe a big speedup, this points at coordinate compatibility checks as the culprit. From there you can experiment with the various options for open_mfdataset, such as coords='minimal', compat='override', etc.

Once you post your file details, we can provide more concrete suggestions.

Wow, thank you!

This is an amazing bug. The defaults say data_vars="all", coords="different" which means always concatenate all data_vars along the concat dimensions (here inferred to be “time”) but only concatenate coords if they differ in the different files.

When decode_cf=False , lat ,lon are data_vars and get concatenated without any checking or reading. When decode_cf=True, lat, lon are promoted to coords, then get checked for equality across all files. The two variables get read sequentially from all files. This is the slowdown you see.

Once again, this is a consequence of bad defaults for concat and open_mfdataset.

I would follow https://docs.xarray.dev/en/stable/user-guide/io.html#reading-multi-file-datasets and use data_vars="minimal", coords="minimal", compat="override" which will only concatenate those variables with the time dimension, and skip any checking for variables that don’t have a time dimension (simply pick the variable from the first file).

OK, so it seems that we need a change to disable wrapping dask arrays with LazilyIndexedArray. Dask arrays are already lazy!

This issue is almost seven years old! It has been “fixed” many times since my original post, but people keep finding new ways to make it reappear. 😆

It seems like having better diagnostics / logging of what is happening under the hood with open_mfdataset is what people really need. Maybe even some sort of utility to pre-scan the files and figure out if they are easily openable or not.

I’ve recently been trying to run open_mfdataset on a large list of large files. When using more than ~100 files the function became so slow that I gave up trying to run it. I then came upon this thread and discovered that if I passed the argument decode_cf=False the function would run in a matter of seconds. Applying decode_cf to the returned dataset after opening then ran in seconds and I ended up with the same dataset following this two step process as I did before. Would it be possible to change one of:

1. where decode_cf is called in open_mfdataset — essentially, open the individual datasets with decode_cf=False and then apply decode_cf to the merged dataset before it is returned;
2. change decode_cf=False to be the default in open_mfdataset?

To me the first solution feels better and I can make a pull request to do this.

From reading this thread I’m under the impression that there’s probably something else going on under the hood that’s causing the slowness of open_mfdataset at present. Obviously it would be best to address this; however, given that the problem was first raised in 2017 and a solution to the underlying problem doesn’t seem to be forthcoming I’d be very pleased to see a “fix” that addresses the symptoms (the slowness) rather than the illness (whatever’s going wrong behind the scenes).

An update on this long-standing issue.

I have learned that open_mfdataset can be blazingly fast if decode_cf=False but extremely slow with decode_cf=True.

As an example, I am loading a POP datataset on cheyenne. Anyone with access can try this example.

base_dir = '/glade/scratch/rpa/'
prefix = 'BRCP85C5CN_ne120_t12_pop62.c13b17.asdphys.001'
code = 'pop.h.nday1.SST'
glob_pattern = os.path.join(base_dir, prefix, '%s.%s.*.nc' % (prefix, code))

def non_time_coords(ds):
    return [v for v in ds.data_vars
            if 'time' not in ds[v].dims]

def drop_non_essential_vars_pop(ds):
    return ds.drop(non_time_coords(ds))   

# this runs almost instantly
ds = xr.open_mfdataset(glob_pattern, decode_times=False, chunks={'time': 1},
                       preprocess=drop_non_essential_vars_pop, decode_cf=False)

And returns this

<xarray.Dataset>
Dimensions:     (d2: 2, nlat: 2400, nlon: 3600, time: 16401, z_t: 62, z_t_150m: 15, z_w: 62, z_w_bot: 62, z_w_top: 62)
Coordinates:
  * z_w_top     (z_w_top) float32 0.0 1000.0 2000.0 3000.0 4000.0 5000.0 ...
  * z_t         (z_t) float32 500.0 1500.0 2500.0 3500.0 4500.0 5500.0 ...
  * z_w         (z_w) float32 0.0 1000.0 2000.0 3000.0 4000.0 5000.0 6000.0 ...
  * z_t_150m    (z_t_150m) float32 500.0 1500.0 2500.0 3500.0 4500.0 5500.0 ...
  * z_w_bot     (z_w_bot) float32 1000.0 2000.0 3000.0 4000.0 5000.0 6000.0 ...
  * time        (time) float64 7.322e+05 7.322e+05 7.322e+05 7.322e+05 ...
Dimensions without coordinates: d2, nlat, nlon
Data variables:
    time_bound  (time, d2) float64 dask.array<shape=(16401, 2), chunksize=(1, 2)>
    SST         (time, nlat, nlon) float32 dask.array<shape=(16401, 2400, 3600), chunksize=(1, 2400, 3600)>
Attributes:
    nsteps_total:  480
    tavg_sum:      64800.0
    title:         BRCP85C5CN_ne120_t12_pop62.c13b17.asdphys.001
    start_time:    This dataset was created on 2016-03-14 at 05:32:30.3
    Conventions:   CF-1.0; http://www.cgd.ucar.edu/cms/eaton/netcdf/CF-curren...
    source:        CCSM POP2, the CCSM Ocean Component
    cell_methods:  cell_methods = time: mean ==> the variable values are aver...
    calendar:      All years have exactly  365 days.
    history:       none
    contents:      Diagnostic and Prognostic Variables
    revision:      $Id: tavg.F90 56176 2013-12-20 18:35:46Z mlevy@ucar.edu $

This is roughly 45 years of daily data, one file per year.

Instead, if I just change decode_cf=True (the default), it takes forever. I can monitor what is happening via the distributed dashboard. It looks like this:

There are more of these open_dataset tasks then there are number of files (45), so I can only presume there are 16401 individual tasks (one for each timestep), which each takes about 1 s in serial.

This is a real failure of lazy decoding. Maybe it can be fixed by #1725, possibly related to #1372.

cc Pangeo folks: @jhamman, @mrocklin

@shoyer I will double check what the bottle neck is and report back. @dcherian interestingly, the parallel option doesn’t seem to do anything when decode_cf=True. From looking at the dask dashboard it seems to load each file sequentially, with each opening being carried out by a different worker but not concurrently. I will see what I can do minimal example wise!

@rabernat wrote:

An update on this long-standing issue.

I have learned that open_mfdataset can be blazingly fast if decode_cf=False but extremely slow with decode_cf=True.

I seem to be experiencing a similar (same?) issue with open_dataset: https://stackoverflow.com/questions/71147712/can-i-force-xarray-open-dataset-to-do-a-lazy-load?stw=2

So is there any word on a best practice, fix, or workaround with the MFDataset performance? Still getting abysmal reading perfomance with a list of NetCDF files that represent sequential times. I want to use MFDataset to chunk multiple time steps into memory at once but its taking 5-10 minutes to construct MFDataset objects and even longer to run .values on it.